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ABSTRACT 


Air-sea exchanges of heat and momentum are important elements in understanding and 
skillfully predicting tropical cyclone intensity, but the magnitude of the corresponding 
wind-speed dependent bulk exchange coefficients is largely unknown at major hurricane 
wind speeds greater than 50 m s'^ Since direct turbulent flux measurements in these 
conditions are extremely difficult, the momentum and enthalpy fluxes were alternatively 
deduced via axisymmetric angular momentum and total energy budgets. A 
comprehensive error analysis was performed using both idealized numerical simulations 
to quantify and mitigate potentially significant uncertainties resulting from unresolved 
budget terms and observational errors. An analysis of six missions from the 2003 
CBLAST field program in major hurricanes Fabian and Isabel was conducted using a 
new variational technique. This analysis indicates a near-surface mean drag coefficient 
(Cd) of 2.4x10'^ with a 46% standard deviation and a mean enthalpy coefficient (Ck) of 
1.0x10' with a 40% standard deviation for wind speeds between 52 and 72 m s' . These 
are the first known estimates of Ck and the ratio of enthalpy to drag coefficient (Ck/Cd) 
in major hurricanes. The results suggest that there is no significant change in the 
magnitude of the bulk exchange coefficients estimated at minimal hurricane wind speeds, 
and the ratio Ck/Cd is likely less than 0.75 for wind speeds greater than 50 m s'\ 
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I. 


INTRODUCTION 


Heat, moisture, and momentum exchange at the air-sea interface are primary 
processes in hurricane intensification and maintenance with important theoretical and 
practical implications (Kleinschmidt 1951; Malkus and Riehl 1960; Ooyama 1969; 
Emanuel 1986, 1995b). Wind-speed dependent bulk aerodynamic formulas often are used 
to represent the turbulent fluxes associated with air-sea interactions, relying on drag (Cd), 
moisture (Ce), and enthalpy (Ck) exchange coefficients. However, as Ooyama (1969) 
stated, “Unfortunately, there is little information on Ce under hurricane conditions, other 
than the semispeculative guess that the exchange coefficients of latent heat, sensible heat 
and momentum are probably of the same magnitude.” Nearly 40 years after that 
pioneering study, enthalpy and momentum exchange coefficients are still largely 
unknown at major hurricane wind speeds (>50 m s'\ equivalent to category three and 
higher on the Saffir-Simpson scale). There are known sensitivities of theoretical and 
numerical models of major hurricanes to the magnitude of these exchange coefficients. It 
is therefore important to improve our understanding and forecasts of these phenomena, 
given the high impact that their damaging winds and storm surge can have on coastal 
populations, global economics, marine, and Naval operations. 

It is very difficult to directly measure hurricane surface layer fluxes over the 
ocean due to the extreme conditions and challenging deployment of instrumentation. In 
situ measurements by ships or manned aircraft are very hazardous, and fixed sensors are 
unlikely to be in the correct location or robust enough to obtain the needed 
measurements. Due to these hazards and difficulties in obtaining direct air-sea flux 
measurements, only fluxes at higher altitudes and minimal hurricane force winds have 
been obtained (Black et al. 2007; French et al. 2007; Drennan et al. 2007; Zhang et al. 
2008). Measurements of the drag and enthalpy exchange coefficients at high wind speeds 
also have been made in laboratory experiments (Donelan 2004; Haus et al. 2010), but are 
limited to wind speeds less than ~50 m s'\ Global positioning system (GPS) 
dropwindsonde profiles were used to calculate flux profiles and deduce the drag 
coefficient at the highest wind speeds yet (Powell et al. 2003; Vickery et al. 2009), but 


1 



this technique has not been attempted for retrieving the enthalpy exchange coefficient. 
Indirect retrievals of the drag coefficient in hurricanes have been conducted using ocean 
measurements (Shay and Jacob 2006; Jarosz et al. 2007), but this methodology is not 
applicable for deducing enthalpy exchange coefficients. 

An indirect approach to retrieving momentum exchange using an angular 
momentum budget was originally proposed over fifty years ago (Palmen and Riehl 1957). 
In this formulation, the surface stress is obtained by measured transports of absolute 
angular momentum in an axisymmetric cylindrical coordinate system. If the flux of this 
quantity is known at the top and sides of a prescribed control volume, the flux at the air- 
sea interface can be obtained via residual. A similar budget can be derived for total 
energy to deduce the enthalpy exchange coefficient. The central focus of this research 
was the application of this budget method to major hurricanes using recent observations. 
After Palmen and Riehl (1957) introduced the budget methodology, additional studies 
used the angular momentum budget to deduce the drag coefficient with improved datasets 
(Miller 1962, 1964; Hawkins and Rubsam 1968; Hawkins and Imbembo 1976). The 
quality and density of hurricane observations has improved significantly since that time, 
as well as the analysis techniques used to composite the data. 

In practical application for this dissertation research, some of the budget terms were 
very difficult to calculate with observational data, and known uncertainties also must be 
acknowledged. To gain insight into the relative importance of the various sources of error in 
the energy/momentum budget method, the sensitivity to errors in unresolved budget terms, 
sea-surface temperature (SST), center placement, gridding method, and size of control 
volume were examined systematically using numerical model data where the surface fluxes 
were known. These results were then utilized to determine a quantitative estimate of the 
confidence level of the magnitude of the retrieved bulk exchange coefficients. A high- 
resolution dataset collected in Hurricanes Fabian and Isabel (2003) as part of the Coupled 
Boundary Layers Air-Sea Transfer (CBLAST) experiment (Black et al. 2007) was used to 
apply this energy/momentum budget method to real tropical cyclones. Both tropical cyclones 
(TCs) achieved category four intensity and the data collected present a significant advance in 
the measurement of major hurricanes. Although comprehensive analyses of the structure of 
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Hurricane Isabel have been performed elsewhere using dropsonde, in situ flight level, 
Doppler radar, and satellite data (Montgomery et al. 2006; Aberson et al. 2006; Bell 2006; 
Bell and Montgomery 2008), one unique aspect of the dataset has yet to be fully utilized. A 
series of rapid dropsonde releases in the eyewall region of these two hurricanes were 
performed to accurately construct control volume composites for use in the budget method. 
These dropsonde “sequences” are used in conjunction with Doppler radar, in situ, and 
radiometer data to construct composite analyses for the air-sea exchange calculations 
presented in this dissertation. 

To deduce momentum and energy fluxes from these sequences, the data must be 
gridded in an axisymmetric coordinate system moving with the tropical cyclone. A simple 
objective analysis scheme (Barnes 1973) was used in the previous Isabel studies to obtain 
the kinematic and thermodynamic structure. Recent advances in data assimilation 
techniques will be used here to improve this analysis by deriving a variational procedure 
that provides a maximum likelihood estimate of the gridded structure given estimates of 
background and observational errors. Additionally, this system allows for the incorporation 
of Doppler radar data, which increases the spatial coverage and improves the kinematic 
measurements in the analysis. This integrated variational technique also can provide a 
unique framework for future theoretical TC studies and numerical model initialization. 

To provide a broad scientific context for the current research. Chapter II is a review 
of relevant previous studies on air-sea fluxes in tropical cyclones. A derivation of the 
equations for total energy and angular momentum conservation in an axisymmetric, 
cylindrical coordinate system necessary for calculations using the budget method is provided 
in Chapter III, and this is followed by a comprehensive error analysis. This error analysis 
includes results from numerical simulations of a hurricane vortex, wherein the sensitivities to 
simulated observational deficiencies are explored. The variational methodology for deducing 
the axisymmetric tropical cyclone stmcture from real observations is described in Chapter 
rv. The results of the analysis using the described methods are presented in Chapter V. The 
concluding chapter summarizes the key findings of this study and discusses potential 
implications for both theoretical understanding of hurricane intensity and numerical weather. 
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II. TROPICAL CYCLONE SURFACE FLUXES 


A. THEORETICAL AND PRACTICAL IMPLICATIONS 


Reynolds decomposition of the momentum and thermodynamic equations leads to 
the well-known turbulence closure problem, in which the number of unknown variables 
exceeds the number of equations for the prediction of atmospheric flow (Stull 1988). In 
the surface layer, where the turbulent fluxes are nearly constant, the fluxes are often 
parameterized by bulk aerodynamic formulae that require wind-speed dependent 
exchange coefficients for momentum, heat, and moisture. In a tropical cyclone, the 
surface stress can be represented in a cylindrical coordinate system in which the 
tangential (swirling) wind is the dominant component as 

^:.e = -p'v'v' = C^p I I V , s 


where v denotes the tangential wind, w the vertical wind, p the density, II the 

horizontal wind speed, Cd the drag coefficient, primes denote perturbation quantities, and 
the overbar is a Reynolds averaging operator. In the axisymmetric coordinate system 
used in this study, the Reynolds averaging operator encompasses both azimuthal and 
temporal averaging. A similar expression for the enthalpy flux can be derived: 

+ ^zc, = + pLw'q' = C^p\ fi, \{k* -k) ^2.2) 

where T denotes the temperature, Cp the specific heat of air, L the latent heat of 
vaporization, q the water vapor, k the enthalpy, k* the saturation enthalpy at the sea 
surface, and Ck the bulk enthalpy exchange coefficient (Emanuel 1995). Since the ocean 
surface characteristics change significantly with increasing wind speed (Black et al. 
1986), a simple extrapolation of the bulk exchange coefficient magnitudes derived at low 
wind speeds to 50 m s‘^ and beyond is not necessarily justified. Emanuel (2003) 
presented a similarity hypothesis for enthalpy and drag coefficients at very high wind 
speeds using dimensional analysis, and proposed that their magnitudes become 
independent of wind speed at approximately hurricane force and higher (>33 m s'^). That 
analysis suggested also that the enthalpy coefficient should vary as a function of 
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temperature, sueh that the ratio Ck/Cd would be approximately unity at typical near¬ 
surface air temperatures in major hurricanes. The earlier work of Ooyama (1969), 
Rosenthal (1971), and Emanuel (1986, 1995) suggests a strong dependence of Ck/Cd on 
intensity. Therefore, validating or rejecting these hypotheses about the behavior of the 
bulk exchange coefficients at major hurricane wind speeds has important implications for 
tropical cyclone intensity theory and numerical weather prediction (NWP). 

The importance of air-sea energy and momentum exchange in hurricanes was first 
introduced by Kleinschmidt (1951) in the form of an energy balance between the amount 
of sensible and latent heat a TC can extract from the ocean surface and the momentum 
lost to the sea through frictional dissipation. Malkus and Riehl (1960) presented a 
potential intensity (PI) theory based on the concept of a steady-state TC at maximum 
intensity that maintains this energy balance, and they derived an expression for the 
maximum tangential wind of a TC given the a priori SST, air temperature, and bulk heat 
and momentum exchange coefficients. While they acknowledged the importance of 
latent heat fluxes from the sea surface, the bulk moisture exchange coefficient was not 
explicitly included in their equation for maximum tangential wind. Ooyama (1969) and 
Rosenthal (1971) were the first to explicitly link the flux of water vapor, or latent heat, 
with numerically simulated hurricane intensity. Emanuel (1986) expanded these concepts 
and derived an expression for maximum intensity assuming a steady-state balance 
between the amount of energy obtained from the ocean and that dissipated by friction in 
the boundary layer, which is similar in general terms to that presented by Kleinschmidt 
(1951) and Malkus and Riehl (1960). The primary parameters that govern Emanuel’s PI 
(E-PI) are the sea-surface temperature, the outflow temperature, boundary layer relative 
humidity, and the ratio of the bulk enthalpy (Ck) and momentum (Cd) exchange 
coefficients (Ck/Cd). Due to limited high-resolution surface flux measurements in the 
extreme conditions of a hurricane eyewall, exchange coefficients are the most uncertain 
of these parameters at wind speeds above category one (Black et al. 2007). 

Bell (2006), Montgomery et al. (2006), and Bell and Montgomery (2008) carried 
out an extensive test of the E-PI predictions using observations collected during 
Hurricane Isabel as part of the 2003 CBEAST field program. Since little guidance existed 
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for the magnitude of the bulk exchange coefficients at major hurricane wind speeds for 
the main potential intensity estimate, the ratio of bulk enthalpy and momentum exchange 
coefficients was assumed to be unity (Ck/Cd =1). The range of E-PI estimates from Bell 
and Montgomery (2008) is ~30 m s'' by changing the assumptions in the calculations 
(Fig. 1, note that all figures appear at the end of the subsections). The upper-bound 
estimate neglects entirely the ocean cooling, includes dissipative heating, and assumes 
that Ck/Cd = 1. The lower-bound estimate assumes Ck/Cd = 0.5 and that the ocean 
cooling effect is compensated by the dissipative heating effect. Accurate determination of 
the bulk exchange coefficients would clearly help narrow the uncertainty of these 
estimates. 

Recent research (Smith et al. 2008; Bryan and Rotunno 2009) has suggested that 
the primary discrepancy between observed and simulated axisymmetric intensities and 
the E-PI theoretical limit and is due to the implicit assumption of gradient wind balance 
in the boundary layer. In three-dimensional simulations by Montgomery et al. (2010), the 
maximum intensity has a much lower sensitivity to the surface drag coefficient than in 
the axisymmetric simulations. The role of radial transport of moist entropy across the 
eye-eyewall interface in two dimensions (Persing and Montgomery 2003) and three 
dimensions (Eastin et al. 2002; Cram et al. 2003) is also not currently included in PI 
theory. While it is clear that a revised PI theory will need to be developed to address 
these legitimate limitations, this does not diminish the importance of estimating surface 
fluxes at high wind speeds as part of accurate determination of hurricane maximum 
intensity. 

The importance of understanding of air-sea exchange at high wind speeds is not 
limited to the theoretical domain. While tropical cyclone track forecasts have been found 
to be relatively insensitive to the parameterization of surface fluxes of enthalpy and 
momentum, the ability to skillfully predict tropical cyclone intensity and structure has 
been found to be dependent on meaningful and accurate parameterization of these 
processes (Davis et al. 2008). Ooyama (1969) was the first to demonstrate the importance 
of parameterized air-sea fluxes for growth and maintenance of a tropical cyclone using a 
simplified three-layer model of the moist atmosphere. Rosenthal (1971) and Rotunno and 
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Emanuel (1987) further demonstrated this importance in hurricane models, and Emanuel 
(1995) demonstrated sensitivity of the simulated intensity to the prescribed surface 
exchange coefficients in a simplified model framework. Subsequent studies by Braun and 
Tao (2000), Bao et al. (2002), and Davis et al. (2008) using full physics models have also 
shown the sensitivity of hurricane intensity forecasts to modifications in the surface 
roughness length parameterization used to calculate heat and momentum fluxes. 
Hurricane simulations using the Navy Coupled Ocean/Atmosphere Mesoscale Prediction 
System are also sensitive to the surface flux and sea spray parameterizations (J. Doyle 
2009, personal communication). Despite improvements in the quality of the intensity 
guidance available to forecasters at the National Hurricane Center, relatively little 
improvement has been achieved in the forecast of tropical cyclone intensity since 1990 
(Rappaport et al. 2009). Reducing errors in the numerical weather model representation 
of air-sea interaction in major tropical cyclones would likely contribute to improved 
intensity forecasts issued by the Joint Typhoon Warning Center and the National 
Hurricane Center. 
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Figure 1. 
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Theoretically predicted azimuthal mean Vmax at the boundary layer top for 
varying outflow temperature and near-core SST with a constant RH = 80% 
(dashed line). ‘X’ indicates the primary potential intensity estimate for the 
observed near environment around Isabel (using 27.5 °C SST near eyewall 
region associated with Fabian’s wake), which yields a 56.6 m s’' mean Vmax- 
The dark solid curve represents the average storm-relative tangential wind 
speed at the top of the boundary layer derived from the dropwindsonde 
measurements. The shading represents the 6 m s'^ standard deviation of this 
mean value. This a priori E-PI estimate assumes Ck/Cd =1. From Bell and 
Montgomery (2008). 
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B. PREVIOUS SURFACE FUUX OBSERVATIONS 


Surface fluxes are typically measured by one of three methods-eddy correlation, 
inertial dissipation, or flux-profile relationships. Each of these methods has different 
advantages and disadvantages depending on the measurement apparatus and field 
conditions, but all three present difficulties in the TCBL. At major hurricane wind speeds, 
the main challenge for these types of measurements is placement and robustness of the 
sensors. Given that the spatial and temporal occurrence of these wind speeds is very 
small compared to the ocean basin size and seasonal timescales in which hurricanes 
occur, the probability of a fixed sensor encountering them is very low. Additionally, 
operational ocean sensors (such as buoys) typically are not designed to withstand the 
extreme conditions found in the open ocean at the base of a major hurricane eyewall, and 
instrumented ships prudently avoid these conditions as well. Therefore, laboratory 
experiments and observations from research aircraft deployed in tropical cyclones are 
currently the most viable choices for obtaining information about surface fluxes at high 
wind speeds. However, low-altitude flying in the turbulent boundary layer at the base of 
an eyewall is also very hazardous, which requires that remote sensing or expendables 
(i.e., dropwindsondes) be used to obtain measurements in this region. This generally 
limits the use of the eddy correlation or inertial dissipation methods except in special 
circumstances where it is safe to fly in the boundary layer and the aircraft is equipped 
with high-resolution sensors. Similarly, the flux-profile method was not viable until the 
development of the modern NCAR GPS dropwindsonde (Hock and Franklin 1999, 
hereafter “dropsonde”). It is perhaps not surprising, then, that the first attempts to 
determine surface fluxes in tropical cyclones were performed by indirect retrievals 
similar to that proposed in this study. 

Palmen and Riehl (1957) reported the first calculation of the drag coefficient in a 
tropical cyclone using a tangential wind budget using two composite wind fields created 
by Hughes (1952) at low levels and Jordan (1952) at upper levels. Hughes (1952) 
analyzed 84 flights from 28 storms in which reconnaissance was available during 1945- 
1947, while Jordan (1952) used rawinsonde observations collected from islands and 
coastal stations near tropical cyclones during the period 1945-1951. These two 
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composites were combined to create an axisymmetric tropical cyclone structure with one 
degree latitude horizontal resolution and variable vertical resolution for use in the budget 
calculations. The original resolution of the low-level data was 0.5 degrees with the lowest 
reported winds at 500 ft (152 m), and no upper-air data were obtained within -200 km of 
the center. To derive the surface stress, the tangential momentum equation was integrated 
from the surface to the top of the inflow layer, assuming steady-state and cylindrical 
axisymmetry, which yielded: 


^--e = 



(2.3) 


where denotes the absolute vorticity / + d{rv)/rdr , and pO and pH are the pressures at 

the surface and top of the inflow layer, respectively. Vertical motion was derived from 
the mass continuity equation with a presumed zero boundary condition at the surface. 
This calculation yielded a Cd of 2.1x10' for the innermost radial ring with a surface layer 
wind speed of -26 m s'\ but it apparent that this was a coarse estimate given the 
extensive averaging used to create the dataset. Subsequent budgets of angular momentum 
and energy were calculated in that study, but no computation of an enthalpy exchange 
coefficient was attempted. 


This budget method was repeated in several subsequent studies with improved 
datasets. Miller (1962) used an aircraft dataset collected in Hurricane Helene (1958) to 
extend the Palmen and Riehl (1957) results by using a single storm as opposed to a 
composite, and at higher surface wind speeds. This study yielded drag coefficients 
ranging from 2.4x10' to 3.2x10' for winds from 30 to 40 m s' . Aircraft data from 
Hurricane Donna (1960) were used to calculate drag coefficients of 3.6x10' and 4.19 
xlO' at 36 and 52 m s' , respectively by Miller (1964), as shown in Figure 2. These 
results indicated that drag coefficient continued to increase with wind speed as storms 
reached major hurricane intensity. Hawkins and Rubsam (1968) and Hawkins and 
Imbembo (1976) continued these analyses with Hurricane Hilda (1964) and Hurricane 
Inez (1966), and extended the wind speed behavior of the drag coefficient to nearly 70 
m s' as shown in Figure 3. To match the maximum drag coefficient of 4.6x10' at 67 

m s'^ found using the Inez dataset, the least-squares curve fit was changed from linear to 
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quadratic. Hawkins and Rubsam (1968) also included a more comprehensive error 
analysis of the budget method, calculating the surface stress using both the integrated 
tangential momentum equation (so-called “vorticity” method) and the angular momentum 
form, which is similar to the current study. These two methods are analytically 
equivalent, but produced slightly different numerical results. Also particularly relevant to 
the current study, Hawkins and Rubsam also compared the results obtained when 
changing the top of the integration volume from the top of the inflow layer to the top of 
the atmosphere (100 hPa). The differences between the two integration methods and 
different boundary conditions suggested an error in the drag coefficient estimates of 
about +30%. An energy budget was also calculated, but the Ch and Ce values were 
assumed to be equal in magnitude to the derived Cd, as opposed to being retrieved from 
the budget itself. In summary, none of these studies using budget methods indicated that 
the drag would level off or decrease at major hurricane wind speed, and none presented 
estimates of a heat or enthalpy exchange coefficient. Frank (1984) also used an angular 
momentum budget to estimate the drag coefficient from a composite analysis of 
Hurricane Frederic (1979). Frank estimated a ±50% error in Cd using this method by 
calculating a standard deviation of estimates from different radial bands. 

Moss and Rosenthal (1975) used the same aircraft dataset from Hurricane Inez 
(1966) to calculate the drag coefficient using the Deardorff boundary layer 
parameterization (Deardorff 1972), and their results compared very favorably to those 
obtained by Hawkins and Imbembo (1976) using the angular momentum budget. They 
also provided an estimate of the bulk heat exchange coefficient at 48 m s'', but they 
admitted that it clearly had a high bias. Moss and Rosenthal also directly compared these 

Q c 

drag coefficients with the simple Deacon’s formula (Cd =1.1x10' +4x10'IV1) and a 
formula proposed by Miller (1969, Cd = 1x10'^ + 7xlO'^IV1), and concluded that these 
equations gave values that were too low at high wind speeds, but too high at low wind 
speeds for Inez, but were reasonable for the dataset from Hurricane Daisy (1958). 
Alhough not specifically conducted in the TCBL, Large and Pond (1981, 1982) 
performed eddy correlation and inertial dissipation estimates of open ocean heat and 
momentum exchange in wind speeds up to ~25 m s'', and their results have been 
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validated by subsequent studies (Smith et al. 1992). These measurements added further 
support to the hypothesis that the drag coefficient continues to increase with wind speed 
up to 25 m s'\ In contrast, Large and Pond (1982) showed no significant wind speed 
dependence on the moisture flux in this wind speed range, and this result has also been 
further validated (DeCosmo et al. 1996). 

Powell et al. (2003) used a flux-profile relationship to infer the drag coefficient 
from dropsonde profiles released in hurricane eyewalls. The wind speed is obtained via 
the displacement of the dropsonde’s GPS position as it descends toward the ocean surface 
after release from a research aircraft. Although the dropsonde is a Lagrangian 
measurement as it travels several kilometers downwind as it falls, the wind profiles can 
be treated in an Eulerian framework by assuming the azimuthal variation in the wind 
speed is small compared to the vertical variation. Fitting the vertical profile to a least- 
squares logarithmic line yields the surface roughness (zq) as the intercept and the friction 
velocity (u*) as the slope from the flux-profile relationship in a neutral surface layer. The 
implied drag coefficient from these results are plotted with estimates from the earlier 
budget studies in Figure 4. Powell et al. provided the first indications of a decrease in the 
drag coefficient from maximum values around -2.5x10' at 30-40 m s' . Additional 
dropsonde profiles have been collected in major hurricanes since the 2003 study, 
including some used in this research, which have extended the wind speed dependence to 
near 60 m s'^ (Figure 5) with a near constant drag coefficient even at this intensity 
(Vickery et al 2009). 

Additional evidence that the drag coefficient does not increase above -35 m s'^ 
was reported by Donelan et al. (2004) using laboratory tank measurements. The eddy 
correlation method was used for wind speeds up to 26 m s'', and a momentum budget 
retrieval was used for wind speeds from 20-50 m s''. These tank measurements showed a 
“saturation” of the drag coefficient around hurricane force wind speed (33 m s''), and 
suggested a limiting aerodynamic roughness of the surface waves above these speeds. 
The momentum budget approach used in the laboratory has also been utilized in field 
experiments with ocean current measurements. Shay and Jacob (2006) used airborne 
expendable current profilers released into Hurricane Gilbert (1988) to estimate the 
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downward kinetic energy flux and retrieve the drag coefficient up to ~40 m s'\ which 
leveled off above 28 m s'^ to a value near 3.5x10'^ (dashed line in Figure 4). A similar 
approach was used by Jarosz et al. (2007) to deduce the drag coefficient from acoustic 
Doppler current profilers moored in the path of Hurricane Ivan (2006). A momentum 
budget using this dataset indicated a quadratic relationship with wind speed, with a peak 
value near 2.5x10'^ at ~33 m s'^ that then decreased to 1.5x10'^ at 45 m s'\ 

The uncertainty in the magnitude of the drag and enthalpy exchange coefficients 
in the TCBL was one of the factors leading to the deployment of the NOAA P-3 aircraft 
as part of the CBLAST field campaign (Black et al. 2007). Due to the safety hazards 
mentioned previously, the aircraft were only flown in the clear air boundary layer 
between rainbands, but fortunately favorable conditions for turbulence observations were 
achieved in two major hurricanes. French et al. (2007) reported the measurements of 
momentum flux from 48 flux calculation legs flown within 400 m of the surface in 
Hurricanes Isabel and Fabian (2003). These results provided the first open-ocean eddy 
correlation measurements in the TCBL at hurricane-force wind speeds. Estimates of the 
drag coefficient versus wind speed from this study (Figure 6) indicated no discemable 
dependence on speed in the range measured. Although the French et al. results slightly 
differ from those reported by Powell et al. (2003) and Donelan et al. (2004), general 
agreement exists that the extrapolated Large and Pond (1981) formula results in too high 
drag coefficients. Latent heat (Drennan et al. 2007) and enthalpy (Zhang et al. 2008) 
fluxes derived from eddy correlation measurements during CBLAST are consistent with 
a lack of dependence of the drag coefficient on wind speed as suggested in other studies. 
Recent laboratory research also indicates that the enthalpy exchange coefficient is nearly 
independent of wind speeds between 13 and 40 m s'^ (Haus et al. 2010, Jeong et al. 
2010). Enthalpy exchange coefficients as a function of wind speed from the HEXOS, 
CBLAST, and laboratory experiments taken from Haus et al. (2010) indicate some scatter 
but no significant trend with increasing wind speed above ~10 m s‘^ (Ligure 7). 

In summary, most recent research has indicated a steady or slightly decreasing 
drag coefficient with increasing wind speed beyond 30 m s"', and also near-constant 
enthalpy exchange coefficient in hurricane-force winds. The limited estimates of the drag 
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coefficient above 50 m s'^ are inconclusive as to the magnitude, with earlier budget 
studies suggesting a substantially higher drag than the flux-profile method of Powell et al. 
(2003) and Vickery et al. (2009). Although no known eddy correlation measurements 
exist at these wind speeds, the laboratory tank momentum budget in Donelan et al. (2004) 
is in closer agreement with the studies of Powell et al. and Vickery et al. The magnitude 
of the enthalpy exchange at major hurricane wind speeds is even more uncertain. While 
there currently is no evidence to suggest that the wind speed independence should not 
continue above 35 m s’', photographs of the sea surface in these conditions depict a 
different character to the sea state than in minimal hurricane winds (Black et al. 1986). 
The sea state at Beaufort category 19 (>50 m s‘^) shown in Figure 8b is described as “low 
clouds, spray, and foam merge into large, white areas frequently referred to as ‘white 
sheets’ by reconnaissance crews” that cover the entire surface at these wind speeds. 
White water only covers about 30-40% of the photograph at ~30 m s‘^ (Figure 8a), with 
narrow, parallel streaks being the defining characteristic. It must be noted that 
parameterizing the complexity of the air-sea interaction at these wind speeds with 10 m 
bulk exchange coefficient may be an over-simplification, but given the established 
theoretical and numerical reliance on these quantities it is still believed to be a useful 
endeavor. 
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Figure 2. 



Drag coefficient Cd as a function of wind speed over water. From Miller 
(1964). 


16 




DRAG COEFFICIENT 


Figure 3. 


DRAG COEFFICIENTS FOR HURRICANE INEZ VS. A BACKGROUND 
OF OTHER HURRICANE AND LOWER SPEED DETERMINATIONS 



WIND SPEED (M/SEC) 


Drag coefficients for Hurricane Inez (1966) from Hawkins and Imbembo 
(1976) plotted with values for Hurricanes Hilda and Helene and some lower 
speed determinations. 
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Variation of the drag coefficient with mean wind speed at 10 m altitude near 
the radius of maximum wind (RMW). From Vickery et al. (2009). 
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Estimated drag coefficient as a function of Uio for the 48 flux runs from 
French et al. (2007) showing the binned values (circles) and the 95% 
confidence interval from this study and extrapolation of results from Large 
and Pond (1981) and Smith (1980), dotted and dashed-dotted, respectively. 
Also shown are results from Donelan et al. (2004; diamonds) and Powell et 
al. (2003; squares) to 42 m s'\ 
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Wind speed dependence of Ck from Haus et al. (2010). ASIST laboratory 
results (•) and CBLAST (A) measurements are shown with HEXOS results 
(x). After binning observations by wind speed, the mean and 95% confidence 
intervals as determined from a t-distribution of the combined HEXOS and 
CBLAST field data are shown in black. 
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Figure 8. Typical appearance of the sea surface in hurricane conditions of Beaufort 
category (a) 11 and (b) 19. The photographs were taken in (a) Hurricane 
Eloise at an altitude of 312 m at 2246 UTC 22 September 1975. The aircraft 
was in the right-front quadrant of the storm about 140 km from the center. 
The flight-level wind was 32.8 m s'^ and the mean and sustained 20 m winds 
were 27 and 30 m s'\ respectively, (b) Hurricane David at an altitude of 454 
m at 1104 UTC 30 August 1979. The aircraft was in the right-rear quadrant 
about 45 km from David’s center. The flight-level wind was 63.0 m s‘\ and 
the mean and sustained winds were 50 and 57 m s'\ respectively. Mean and 
sustained winds are ten and one-minute averages, respectively. From Black et 
al. (1986). 
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III. BUDGET METHODOLOGY 


A. CONSERVATION OF ANGULAR MOMENTUM 


To determine the bulk momentum exchange coefficient, the flux form of the 
azimuthally-averaged tangential momentum equation in cylindrical coordinates is: 


3(pv) , dipruv) ^ dipwv) , puv ^ ^ 

-+-:-+ —:-+- + pfu = Fg 


dt 


rdr 


dz 


(3.1) 


where u is the radial wind, v the tangential wind, w the vertical wind, p the density,/the 
Coriolis parameter, and Eg molecular friction/diffusion. Multiplying Equation (3.1) by the 
radius r, and assuming that Fg is small above the viscous sublayer yields: 


dt 


(prv)+ 


— (pMrv)+ puv 
dr 


+ —{pwrv)+ {pruf}= 0 
dz 


(3.2) 


By the chain rule, the term in brackets reduces to: 


dipujrv)) 

dr 


+ puv ■ 


d{rpu{rv)) 

rdr 


(3.3) 


Considering the final term in braces, the axisymmetric mass continuity equation is: 

dp ^ d(pru) ^ djpw) ^ ^ 

dt rdr dz (3.4) 

Assuming an /-plane approximation, the Coriolis parameter is a constant and therefore 
multiplying the continuity equation by r f 12 yields: 

r^^(Pf) ^ rd(pruf) ^ r^d(pwf) ^ ^ 

2dt 2dr 2dz H 31 


By the chain rule, the second term can be expressed as: 


rd(pruf) 

2dr 


djprufir^)) 

2rdr 


-pruf 


(3.6) 
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Therefore: 


pruf = ^ ^ipwfr^) 

2dt 2rdr 2dz 


(3.7) 


Substituting (3.3) and (3.7) into (3.2) yields: 

|[p(rv + > 72 )] + ^[pru (rv + > 72 )] + jJ_Pw (rv + >72)] = 

d(pM) ^ d(pruM) ^ d(pwM) ^ 
dt rdr dz 


where the absolute angular momentum is defined as M = rv +— fr 


2 


Each variable is then separated into a mean and perturbation quantity, such that 
V = V + v' ^ where the mean is both an azimuthal and temporal average. Reynolds 
averaging and making the Boussinesq approximation (e.g., neglecting p') yields: 


d(pM) d(pruM) d(pwM) d(rpu'M') d(pw'M') 
dt rdr dz rdr dz 


The Reynolds stresses are then defined as = -pu'v' and = -pw'v', where the 
subscripts represent the wind directions of the covariances (r is the radial direction, 6 is 
the tangential direction, and z is the vertical direction). Expanding the last two terms 

d(rpu'M') d(r^pu'v') 3(r>^g) 

rdr rdr rdr 

d{pw'M') d{rpw'v') rdr^g 

dz dz dz (3.10) 


yields the axisymmetric angular momentum equation in flux form: 

d(pM) ^ d(pruM) ^ d(pwM) djrhj 
dt rdr dz rdr dz 


(3.11) 
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Integrating over a control volume from zi to zi and ri to r 2 yields the integrated flux form 
of the conservation of angular momentum in axisymmetric cylindrical coordinates: 


n 

•fz, 


■rdrdz 


+ 


d(pM ) 

^1 *’'^1 dt 

p djpruM - r^z [ 3(pwM - ) 

Jzi Jfi Jr, J; 


dr 


dz 


dz rdr ■■ 
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(3.12) 


Rearranging (3.12) and dividing by 2 k, the radially-integrated turbulent 
momentum flux at zi is then: 

-r^ \puM - rr^g ] dz + r, [puM - ri^g ] dz 
- J ^ [ pwM - rdr -I- 1 ^ [pwM ] 1^^ rdr “ | " J " — ^-rdrdz 


dt 


(3.13) 


From (3.13), the integrated surface stress for an arbitrary control volume may be 
calculated. A schematic of a hypothetical control volume is shown in Figure 9, which 
illustrates the approximate location of zi, Zz, rj, and rz for an idealized hurricane flow. 
Note that zi does not necessarily have to be at 10 m altitude, but could be anywhere in the 
surface layer assuming the fluxes are nearly constant in that layer. 

Assuming the vertical stress at zi can be represented by a bulk 

formula, ~ ^ ~ ^r>P ^ I v ^ neglecting the radial variation of Cd over the 

control volume yields the drag coefficient: 
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Cd= (-'2 ] I, dz + J [ puM - rz^, ] I,, dz 

-j'‘ [p«-M - , r^.] 1.^ rdr + J" [pn M ] I. rcl, - J" J; ^^^rdrdz] / 

J'’[pl«,,li-]l rdr 


(3.14) 


Since the tendency and additional flux terms cannot be calculated using the 
current dataset, the working form for use with these observations is 


Cd = (-'2 [P»3^ ] I, dz + r, j'l' [puM] dz 

- j'' [pvvM ] rdr + J'" [ pu’M ] I. rdr] / (J'' [p II v] L, r\/r) + R 


(3.15) 


where the unresolved residual is 


R = |r 2 [ rz,., ] dz - 1 ] [ rz,., ] ^ dz + [ rz.g ] rdr - ^^^p-rdrdz 

j'' [p I u, I v] rV/r 


(3.16) 


Alternatively, one can obtain the average surface stress by dividing (3.13) by the 

pr2 2 

integrated square of the radius r dr, or divide by the mass-weighted term 

•''■1 

r'2 2 

pr dr to obtain the average friction velocity m* in the control volume. The residual 

•''•1 

terms were also unresolved in the previous studies utilizing a budget methodology. An 
accurate determination of the surface stress, friction velocity, or drag coefficient from this 
formulation requires minimal errors in the specification of the axisymmetric mass and 
wind fields, strict adherence to axisymmetric mass continuity, and a small magnitude of 
the unresolved residual terms. Each of these error sources is examined in Chapter III.D. 
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Figure 9. Schematic illustrating hypothetical control volume (black dashed line) used 

for the budget methodology. A simplified secondary circulation (gray 
streamlines) and region of maximum wind (Vmax) are shown to indicate the 
control volume encompasses the eyewall region. 


B. 


CONSERVATION OF ENERGY 


To determine the bulk moist enthalpy exchange coefficient, the first law of 
thermodynamics in material form with enthalpy as the state variable is: 


T>(C/) ^ ^ 

Dt Dt * ^ ^ 


(3.17) 


where Cp is the specific heat at constant pressure, T is the temperature, k is the thermal 
conductivity, Qr is the diabatic heat exchange due to radiative transfer, Qf is the heat 
arising from frictional dissipation of kinetic energy, and Ql is the latent heat release 
associated with phase changes of water, with the additional terms defined above in 
Chapter III.A. The latent heat release is given by: 


Ql = 


-pL Dcu 
Dt 




(3.18) 


with L the latent heat of vaporization, q is the water vapor mixing ratio, and qsat is 
the saturation mixing ratio. Equation (3.18) can be approximated by 
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Qi^ ~-LpDq / Dt since qsat is small compared to unity, and Ql = 0 when Dq/Dt = 0 
(Gill 1982). The heat arising from frictional dissipation is: 





dx 


V ) 


(3.19) 


where v is the kinematic viscosity, and the indices i and j represent summation over the 
three spatial dimensions. Equation (3.17) then becomes: 


DT ^ Dq Dp 

pc^ + pL — = — + pv 

Dt Dt Dt 




+ kW^T + Q, 


(3.20) 


where the mass continuity equation has been used to put the first two terms into flux 
form. This equation describes the internal (sensible and latent heat) energy of an air 
parcel. To derive the total energy equation, expressions for mechanical (kinetic and 
potential) energy are also needed for the parcel. Starting with the momentum equation in 
vector form: 


^ —+ 2Qxh 
Dt 


= -Vp - pVO + 


(3.21) 


where u is the three-dimensional wind vector, 2Q x « is the Coriolis acceleration, g is 
gravity, and is molecular friction/dissipation. Dotting u into (3.21) yields the kinetic 
energy equation, where Fu has been separated into its diffusive and dissipative parts: 
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(3.22) 


with p denoting the molecular viscosity. The potential energy is given by the 
geopotential: 


P -= PM’VO = wgp 

Dt 


(3.23) 
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Since O depends only on z, the partial derivative with respect to time is zero. Adding 
(3.23) to (3.22) yields the mechanical energy equation: 


1 D{ir + v- + w-) 

-p— - - 

2 Dt 
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(3.24) 


where the pressure advection has been replaced using the definition of the material 
derivative, DpfDt = dpldt + u-Vp. Adding the internal energy equation (3.20) and the 
mechanical energy equation (3.24) yields the total energy equation in material form: 
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(3.25) 


Note that the dissipative heating terms have cancelled since they represent a conversion 
between mechanical and heat energy. Here, it is assumed that the diffusive component of 
friction, thermal conductivity, and radiation on the right side of (3.25) are small and can 
be neglected in this application. The time derivative of p on the right-hand side does not 
cancel because the thermodynamic equation was expressed in terms of the enthalpy 
instead of the internal energy (e.g., c^T = c^T + p/p). In this application, local changes 

in pressure are associated with acoustic waves in the atmosphere and this term is small 
and also neglected. With the right side set to zero. Equation (3.25) is known as 
Bernoulli’s equation (Gill 1982). Expanding the material derivative to flux form in 
azimuthally-averaged cylindrical coordinates: 
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(3.26) 


Each variable is then separated into mean and perturbation quantities as before, 
where the mean is both an azimuthal and temporal average, and are Reynolds-averaged 
using the Boussinesq approximation. For simplification, several terms may be grouped 
and defined: 


.2 , ,.2 , ,,, 2 ^ 



(3.27a) 

-p{u''^ +v'^ +w''^) 


2 

(3.27b) 

F^t= pc pU'T' 

(3.27c) 

Fzt= pc pW'T' 

(3.27d) 

Kq = PLii'q' 

(3.27e) 

K, = 

(3.27f) 


where E is the total energy of the mean flow, e is the turbulent kinetic energy, and the F 
terms are radial and vertical fluxes of temperature and moisture in which the subscripts 
indicate the covariances. The Reynolds averaging yields: 
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(3.28) 


Note that terms involving T-jdu. / dxj cancel due to conversion of mean kinetic energy to 

turbulent kinetic energy (e.g., see Lindzen 1990, 90-92), but flux gradients involving 
3(TyM, ) / dxj remain. For lack of an established term, this is coined the “shear flux,” as it 

represents a loss of kinetic energy through the interaction of turbulent momentum fluxes 
and shearing flow. Grouping terms and integrating over the control volume: 
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(3.29) 


Rearranging (3.29) and dividing by 2n, the radially-integrated enthalpy flux at zi 


is: 
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(3.30) 


It is then assumed that the vertical moist enthalpy flux at zi can be represented by a bulk 
formula: 


F^. = F, + F, ^ pc^yr + pLw'q' = C,p I u, \{k*-k) 


(3.31) 


where k* is the saturation moist enthalpy at the sea surface (k* = + Lq*), and < 7 * is 

the saturation mixing ratio at the surface. Neglecting the radial variation of Ck over the 
control volume yields the bulk enthalpy exchange coefficient: 
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(3.32) 


Similar assumptions as in the momentum equation are made to group unresolved terms 
into a residual: 
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Ck = [p^E][, dz - r, [puE\ dz+V [pw£]l,^ rdr-\ " [pwE}\,^ rdr + E^jj^ar / 
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(3.33) 


where the transport of energy by the interaction of turbulence and the mean wind through 
the lower boundary = J +^'^^ 0 ] I ^dr is written separately since it may be 

Jrj 1 

resolved indirectly via the results from the drag coefficient calculation. Alternatively, one 
can obtain the average enthalpy flux by dividing (3.30) by the integrated radius. The 
unresolved residual is 
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(3.34) 


As with the momentum budget, accuracy in the retrieval of the enthalpy exchange 
coefficient requires: (i) minima! errors in the specification of the axisymmetric 
thermodynamic and kinematic fields; (ii) strict adherence to mass continuity; and (iii) a 
small magnitude of the unresolved residual terms. The validity of these assumptions is 
presented in the following section, using numerically simulated tropical cyclones’ 
budgets. 

C. PROOF OF CONCEPT: VERIFICATION USING SIMULATED DATA 

To test the accuracy of the methodology described in the previous chapter, 
comprehensive momentum and enthalpy budget analyses using idealized numerical 
simulations were undertaken. The first objective is to retrieve bulk drag and enthalpy 
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exchange coefficients with known magnitudes using numerically simulated observations 
of tropical cyclones. The second objective (to be described in Chapter II.D) is to use the 
simulated data sets in a budget error analysis. 

Two numerical weather models were used to test the budget methodology. The 
first was the two-dimensional, axisymmetric hurricane model developed by Rotunno and 
Emanuel (1987, hereafter RE87); and the second was the three-dimensional Advanced 
Research Weather and Eorecasting Model version 3.0.1.1 (hereafter, WRE). Since the 
main goal of the numerical modeling effort was to produce idealized, but still physically 
realistic, hurricane simulations, an effort was made to closely replicate the initial 
conditions and parameterizations of the two models. The initial atmospheric conditions 
used for both models was the thermodynamic background defined by the Jordan (1958) 
mean tropical Atlantic sounding that was modified to be in thermal wind balance (Smith 
2006) with the initial tangential winds specified from an analytic mesoscale vortex 
(RE87, Equation 37). The sea-surface temperature was fixed in the RE87 model at 26.5° 
C, and at 28.0° C in the WRE simulation. Both simulations used warm rain microphysics, 
with a single category for all liquid water but variable terminal fall speed depending on 
the mixing ratio in RE87, and separate categories for cloud and rain water in the WRE 
model using the Kessler (1969) microphysics scheme. The RE87 model used Newtonian 
cooling, and radiation was not included in the WRE simulation. In RE87, a uniform 1 km 
horizontal (radial) and 250 m vertical grid spacing were used in a 1500 km by 25 km 
domain. All other settings are similar to those from the 4X run of Persing and 
Montgomery (2003), except the surface fluxes as described below. In WRE, a four-way 
nested domain was used with a fine mesh containing 205 x 208 gridpoints at 1 km 
horizontal grid spacing, which tripled with each successive mesh to a 5400 km square 
coarse domain with 27 km resolution. Eifty vertical levels were unevenly distributed with 
the highest resolution in the boundary layer using the default WRE sigma levels. The 
WRE model was run for a total of 8 days, and reached steady state around 150 hours into 
the simulation. The RE87 model simulations reached steady state from 100-150 hours 
and were extended to 12 days due to the lower computational requirements. A summary 
of the model settings is given in Table 1. 
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In addition to the fundamental geometrical and dynamical differences of the two 
models, the surface flux and sub-grid turbulence parameterizations use different 
approaches, which provided independent storm structures and magnitudes of the surface 
exchange for testing the budget methodology. In the RE87 model, sub-grid turbulence is 
parameterized by a local, first-order closure that relates turbulent fluxes to resolved grid- 
scale gradients by an eddy-viscosity assumption. Since the model is axisymmetric, all 
three-dimensional motions are therefore parameterized by gradients in the azimuthal 
mean flow. The resolved deformation flow and a prescribed mixing length determine the 
eddy viscosity following Smagorinsky (1963). 

Surface fluxes in the RE87 model are given by bulk aerodynamic formula, such 
that the bulk surface exchange coefficients can be prescribed exactly. In an effort to 
determine whether the proposed budget methodology can distinguish effectively between 
different magnitudes of the drag and enthalpy exchange coefficients, five simulations 
with different coefficient magnitudes and different Ck/Cd ratios were performed. The 
“high” magnitude simulations used the Deacon formula, as originally prescribed in RE87: 

C „ = 1.1x10-V4x10--’(//- + v -)"-|^/2 

with the corresponding Ck given by a defined ratio of either 1.0 (“highl”) or 0.7 
(“highO.7”) to the drag coefficient. The “mid” magnitude simulations used the same 
Deacon formula at surface wind speeds up to 35 m s'^ but the exchange coefficients then 
were capped at higher wind speeds. Note that the surface fluxes were not capped above 
hurricane force, merely the wind-speed dependence of the exchange coefficients. Two 
“mid” simulations using Ck/Cd of 1.0 (“midi”) and 0.7 (“midO.7”) were performed. A 
single “low” simulation also was performed in which the bulk exchange coefficients were 
derived from direct aircraft measurements of the turbulent fluxes reported in Black et al. 
(2004), Drennan et al. (2004), and Zhang et al. (2008) (see Chapter II.B for a review of 
these results). Eor the “low” simulation, Ck was set to 1.1x10' for all wind speeds, and 
the Deacon formula Cd was capped at 1.9x10'^. These five RE87 simulations give a range 
of coefficients for testing the sensitivity of the budget methodology. 
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The WRF model has several turbulence closures available. The modified Yonsei 
University Scheme (hereafter YSU) (Noh et al. 2003; Hong et al. 2006) was selected for 
this study. This scheme employs a non-local boundary layer parameterization in which a 
parabolic eddy viscosity profile is prescribed to match the estimated model boundary 
layer depth (defined in the model as the top of a well-mixed layer of virtual potential 
temperature). Additional vertical fluxes that account for transport by large eddies and 
entrainment at the top of the PEL are also included in this formulation. Above the 
boundary layer, the YSU scheme uses a local eddy viscosity assumption which is 
somewhat similar to the RE87 model except it is formulated in three dimensions. Davis et 
al. (2008) described modifications to the original YSU scheme for hurricane simulations 
in a similar manner to the “mid” RE87 runs described above, such that the drag and 
enthalpy exchange coefficients do not increase beyond 35 m s'\ In Davis et al. (2008), 
the bulk exchange coefficients were not modified directly, but rather were adjusted 
through the surface roughness by setting z,, = 10exp(-10 /m!^‘^) with lower and upper 

limits of 1.25 x 10' m and 2.85 x 10' m, respectively. This modified surface flux option 
was selected for this study, but no modifications were made to the WRE source code. The 
drag and enthalpy exchange coefficients were then inferred from the output u* and 
sensible and latent heat fluxes at each gridpoint. The retrieved WRE fluxes were slightly 
different from those described in Davis et al. (2008), because examination of the v3.0.1.1 
code shows a modified roughness formulation that saturates at 28 m s'^ instead of 35 
m s'^ (the -10 in the exponent is equal to -9 in this version). Although the WRE model 
includes an optional one-dimensional ocean mixed layer model for these simulations, it 
was not included in this study. 

The maximum low-level winds for the simulations are shown in Eigure 10. The 
RE87 winds are reported at 125 m altitude, and the WRE winds are from a diagnostic 10 
m value. All six simulations have instantaneous low-level wind speeds exceeding 70 
m s'\ so that all of these modeled storms are at category five intensity. After a ~4-5 day 
intensification period, all the simulations reach a relatively steady maximum intensity 
ranging from 75—120 m s'\ The weakest storm was the “highO.7” simulation, and the 
strongest was the “midi.” The three-dimensional WRE simulation and the “low” RE87 
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simulation were the two slowest developers, with both storms peaking at 150 hours near 
100 m s'^ before a slight weakening period. The “0.7” Ck/Cd ratio simulations are weaker 
than the corresponding “1.0” simulations for the same flux formulation, which is in 
accord with the predictions of E-PI theory. However, the large difference between the 
“midi” and “highO.7” simulated intensity is counter-intuitive in terms of E-PI theory, 
which depends only on the Ck/Cd ratio and not on the absolute magnitudes of the 
exchange coefficients. The RMW is ~5 km smaller in the “midi” simulation than in the 
“highl” (not shown), which is consistent with a higher tangential wind for a similar 
absolute angular momentum distribution, but a more detailed comparison of the 
differences in the model simulations is beyond the scope of this study. While some of the 
simulated storms are significantly stronger than observed tropical cyclones, they are 
sufficiently realistic for testing the momentum and energy budget methodology outlined 
in Chapters II.B and II.C, with the caveat that the exact magnitudes of various terms may 
be over- or under-estimated compared with actual storms. 

A representative sample of the simulated structures from the RE87 simulations re¬ 
analyzed by the SAMURAI software developed for this study is shown in Eigure 11. The 
SAMURAI analysis is described in detail in Chapter IV. The analysis is shown here to 
facilitate a better comparison of the structure with the observed storms presented in 
Chapter V. A comparison of the raw model results with the SAMURAI analysis is 
discussed in Chapter II.D.4 below. The 6-hour averaged axisymmetric tangential wind 
(Eig. 11, top) at 168 hours from the “low” simulation has a maximum of -110 m s'^ at 
1 km altitude and 20 km radius, with a well-defined secondary circulation. The top of the 
inflow layer has a steep downward slope inside the RMW that approximately matches the 
level of the maximum tangential wind, and outside the RMW the top is at a relatively 
high altitude of 2 km. Upward motion and outflow maxima are found just above the 
maximum tangential wind. The secondary circulation is aligned closely with the angular 
momentum and energy contours in the eyewall region, with a relative energy maximum 
(-351 kJ kg'^) near the surface just radially inward of the RMW. The model vortex has a 
similar kinematic and thermodynamic structure as in previous RE87 simulations (Persing 
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and Montgomery 2003). Although some variations in the structure are simulated 
depending on the flux configuration, all runs are qualitatively similar. 

To obtain the WRF axisymmetric structure, a SAMURAI analysis was performed 
using circulation centers derived using the method described in Chapter V.C.2. Since the 
simulation was on an /-plane with no environmental flow, the wind, pressure, and 
circulation centers were nearly coincident and practically stationary over the period. 
Hourly analyses were created using the 15-minute model output with independent 
estimates of the center position from alternating east-west and north-south slices 
analogous to the “figure-4” flight pattern used in aircraft reconnaissance. No noise or 
center position perturbations were added, which thus yields the “best” analysis possible 
for these simulated fields. The intent was to emulate the analysis procedure used with the 
CBLAST dataset and provide a best-case scenario for understanding errors associated 
with neglecting the residual budget terms discussed in Chapters III.A and III.B, as well as 
small errors introduced in the analysis process. 

An example of the model output used in the analysis is shown in Figure 12, with 
the 1 km altitude fine-mesh horizontal wind vectors and speed and simulated 50 dBZ 
radar reflectivity contour at 144 hours. The eyewall is apparent as a nearly symmetric 
15-km wide annulus of high reflectivity co-located with a region of strong winds over 
80 m s‘^ and a peak wind speed exceeding 120 m s‘^ on the inner edge of the high-rain 
region. Some radial outflow is also evident in the wind vectors at this inner reflectivity 
edge. Minimal structural changes are simulated in the eyewall region over the 24-hour 
period (144-168 hours) used for the budget analysis. 

The axisymmetric SAMURAI analysis at 144 hours for the WRF simulated 
structure (Figure 13) is generally similar to the RE87 structure, but is notably different in 
several details. For example, the wind maximum in the WRF simulation is found at lower 
altitudes and radially inward of the location in the RE87 “low” simulation. Eurthermore, 
the sloping inflow reaches only 1-km height in the outer part of the domain, and is ~20 
m s'^ stronger than the RE87 inflow near the surface, and is topped by a stronger outflow 
region. Mostly congruent secondary circulation and angular momentum and energy 

contours are also found in this simulation, but the energy values are higher. It is noted 
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that the density is lower in the WRF simulation than in the RE87 simulation (not shown), 
such that mass-weighted energy is not as different as Figures 13 and 14 might suggest. 
However, the axisymmetric kinematic and thermodynamic structures of the RE87 and 
WRF simulations are clearly distinct from a budget standpoint. While this is an 
interesting research result, the objective of this section is to test the budget method with 
reasonably realistic, balanced terms, and not to establish the fidelity of the simulations to 
a real tropical cyclone. A systematic exploration of the differences between the 
axisymmetric and 3D models is therefore deferred to a later study. 

Bulk exchange coefficients diagnosed from the model output for the “highl,” 
“midi,” and “low” RE87 experiments and WRF are shown in Figure 14. The RE87 code 
was modified to output all of the subgrid terms, including the residual terms and surface 
stress, heat, and moisture fluxes. Since these surface flux terms were prescribed by the 
bulk aerodynamic formula, the diagnosed exchange coefficients match the analytic values 
very closely, as expected. For the WRF simulations, the equivalent drag and enthalpy 
exchange coefficients were diagnosed using the bulk aerodynamic formula and the output 
friction velocity and surface heat fluxes. These Ck and Cd values are slightly lower than 
the 10 m analytic values since the fluxes were calculated at the lowest sigma level (~40 m 
altitude). For the WRF, the drag coefficient (Fig. 14a) levels off at 28 m s'\ with an 
increasing Ck (Fig. 14b) near 1.9x10' for wind speeds around 60 m s' . This flux profile 
corresponds reasonably well to the “midi” flux configuration used in the RE87 
experiments. All of these magnitudes are within the range of expected values for the 
exchange coefficients in real tropical cyclone boundary layers. 

Since all of the subgrid scale terms could be recorded from the RE87 simulations, 
it was also possible to test the proposed retrieval method when all terms in the budget 
were included. This test for a representative sample from the “highl,” “midi,” and “low” 
simulations is shown in Figure 15. The exact match of the derived and prescribed 
coefficients is as expected, and serves as a validation of the retrieval method when all the 
terms are known. These results were independent of the control volume used for the 
retrieval. Since the goal was to simply validate the retrieval methodology, an equivalent 
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test was not performed with the WRF model due to the added eomplexity of ealeulating 
and storing the Reynolds averaged fluxes in eylindrieal eoordinates. 


Table 1. Summary of numeiioal simulation eonfigurations 


Model 

RE87 

WRF 

Horizontal grid spaeing 

1 km radial 

1 km fine to 27 km eoarse 

Vertieal grid spaeing 

250 m 

50 uneven sigma levels 

starting at 40 m spaeing 

Mierophysies 

Warm rain 

Kessler 

Surfaee Layer 

“highl“highO.7,” “midi,” 

YSU similarity with 


“midO.7,” “low” bulk 

hurrieane modifieation 


aerodynamie formulations 

(isfteflx=l) 

PEL 

L‘ order, loeal elosure 

YSU seheme 

Lateral Boundary 

Jordan sounding 

Jordan sounding 

eonditions 



SST (C) 

26.5 

28.0 

Radiation 

Newtonian eooling 

None 
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Figure 10. 



Simulated maximum low-level wind for the control runs used in sensitivity 
testing. Low-level winds are maximum instantaneous 10 m winds diagnosed 
every 15 minutes from the WRF simulation, and are maximum instantaneous 
axisymmetric 125 m winds every six hours from the five RE87 simulations 
(see inset). 
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RE87'low' 168 hour tangential wind (color), radial wind (contour), and secondary circulation (vector) 



< ^ 


Radius (km) 


8 16 24 32 40 48 56 64 72 80 


0.25 0.5 0.75 1 


96 104 112 120 m S'^ 


50 

10 m/s -► 


1.25 1.5 1.75 2 2.25 2.5 2.75 3 3.25 3.5 3.75 10^fTl^S'^ 


Figure 11. RE87 “low” simulation axisymmetric structure at 144 hours. Top panel 
shows tangential wind (color, m s‘^), radial wind (5 m s'^ contours), and 
secondary circulation. Bottom panel shows angular momentum (color, 
lO^m^s’'), total energy (350 kJ kg'' + 1 kJ kg'' contours), and secondary 
circulation (vector). Solid contours indicate positive values, and dashed 
contours indicate negative values. 
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Figure 12. Wind vectors (scale vector at bottom right) and speed (color scale, m s"') at 1 
km altitude and simulated radar reflectivity 50 dBZ contours at 144 hours 
from the WRF simulation. 
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WRF 144 hour tangential wind (color), radial wind (contour), and secondary circulation (vector) 
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WRF axisymmetric structure at 144 hours. Contours as in Figure 11. 
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Figure 14. 


Diagnosed (a) Cd and (b) Ck from five RE87 simulations and WRF 
simulation (see insets for color symbols) used in the sensitivity analysis. 
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Figure 15. Retrieved exchange coefficients from selected RE87 simulations (see inset) 
using budget methodology with all terms included. 


D. ERROR ANALYSIS 

1. Impact of Unresolved Budget Terms 

The baseline validation of the methodology indicates excellent agreement 
between the derived and prescribed exchange coefficients when all budget terms were 
included. However, several of the terms cannot be estimated with the CBLAST dataset. 
Although the actual magnitude of the unresolved terms in nature is unknown, first-order 
estimates obtained from the RE87 output are shown in Figure 16. The root mean square 
values were determined from 1008 samples of the 6-hourly averaged simulated values 
from 168-288 hours in the RE87 integrations with different flux configurations. Control 
volumes were varied in size and shape from 18.5-44.5 km radius, and from the lowest 
model level (125 m) to 1.625-3.125 km altitude to get a representative sample of the 
various quantities. This analysis indicates that the leading term in the budget is the flux of 
mean quantities through the outer surface of the control volume, and the second leading 
term is the mean flux through the top of the volume—both of these terms are known 

quantities. Although the subtraction of several large-magnitude terms to obtain a small 

45 




residual generally is problematic in finite precision mathematics, in this case the desired 
surface fluxes are the third leading term, so there is some confidence in the ability to 
retrieve this quantity with reasonable accuracy. Radial eddy fluxes were the smallest 
terms in both budgets, which suggests that they can be neglected with minimal error in 
the inner core region. The two mean fluxes across the inner and bottom interface are 
known, with the latter being a smaller term that mainly satisfies integrated mass 
continuity. Unresolved terms of numerical significance are the volume-integrated 
tendency/storage term, the “shear flux” term, and vertical eddy fluxes. Each of these is 
discussed in more detail below. 

Whereas the simulated storms were relatively steady in their intensities over 
multi-day timescales, the wind speed timeseries (Figure 10) indicates that there was non- 
negligible variability on hourly or more frequent timescales. Since the tendency term is 
integrated over the entire control volume, it has the potential to be a significant term in 
the budget. To examine the variability of the tendency term more carefully, budget 
retrievals were performed using both instantaneous and time-averaged model output. A 
finite-difference calculation for the momentum and energy tendencies on the original, 
staggered grid from consecutive model timesteps was found to be close to the 
recalculated value from the momentum (3.11) and energy (3.28) equations using 
unstaggered model output. A power spectrum of the intensity time series has a spectral 
peak at the 25-minute period (not shown), which is consistent with the period of a local 
inertial oscillation (Shapiro and Franklin 1995). However, the magnitude of the 
oscillation may be larger in the axisymmetric simulation than in observed storms. 
Experiments with different averaging periods appeared to damp this oscillation, but 
introduced additional errors in the storage term. Most notably, the unstaggered, time- 
averaged output did not explicitly enforce mass continuity with the numerical precision 
of the staggered, instantaneous gridpoint values. This discrepancy proved to be a 
significant source of error, with fluctuations in the retrieved exchange coefficient 
magnitudes exceeding 500%. This unacceptable level of uncertainty could be reduced by 
strictly enforcing mass continuity through a recalculation of vertical velocity from the 
given divergence field, but this adjustment led to unbalanced fields that required non- 
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negligible tendency terms to close the budget. In this sense, a larger tendency term can be 
alternately viewed as an artificial, residual storage of either momentum or energy inside 
the volume resulting from errors in the resolved fluxes across the top and sides of the 
control volume. These volume-integrated tendency terms also were found to be spatially 
coherent, and were manifest primarily as a sensitivity to the size and shape of the control 
volume if the terms were neglected. A more detailed examination of the sensitivity to the 
control volume is presented in Chapter III.E, but this order of magnitude analysis 
indicates that the estimated error in the exchange coefficients from neglecting this term is 
±20-100% for moderate Cd and Ck values (e.g., between 1 and 2.5 xlO'^). 

The “shear flux” term represents an integrated loss of kinetic energy through the 
lower boundary resulting from the interaction of the mean wind shear and turbulent 

stress, = J ^ + vT^g~\ 1^ rdr. The shear flux term is -40-60% of the magnitude 

of the surface fluxes, which results in a significant low bias if neglected. However, the 
surface stress can be estimated by the bulk aerodynamic formula using the drag 
coefficient derived from the momentum budget, and the mean wind. This correction 
works well with the numerically modeled budgets when the drag coefficient is specified 
exactly, but in the real data the errors due to uncertainties in the surface stress and mean 
wind causes uncertainty. Sensitivity tests indicate that a ±50% error in Cd translates to a 
±20% error in the magnitude of Ck from this term. To avoid adding too much noise to the 
Ck estimates, the mean Cd derived from all six missions was used in this study to 
estimate the surface stress used in evaluating the “shear flux” for Ck- Including errors in 
the mean wind, ±20% appears to be a reasonable estimate of the Ck error introduced by 
estimating this term. 

Unresolved vertical turbulent fluxes at the top of the control volume are typically 
the same sign as the surface fluxes in the lower troposphere, which represents an 
unresolved flux of momentum into the volume, and a flux of energy out of the volume. 
This leads to an underestimation of both the drag and enthalpy exchange coefficients if 
neglected (see Equations 3.14 and 3.32). In the RE87 simulations, these eddy terms 
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showed a strong dependence on the altitude of the volume top, and were 10-30% of the 
surface flux magnitudes for the 1.625-3.125 km tops used in Figure 15. 

For the CBLAST analysis, the choice of volume top is not straightforward. 
Traditional formalism for the planetary boundary layer assumes that the magnitude of the 
turbulent fluxes decreases to zero at the top of the PBL, but an examination of the vertical 
structure of turbulence using the CBLAST dataset (Zhang et al. 2009) suggests a 
separation of the TC boundary layer height defined mechanically and thermodynamically 
(Smith and Montgomery 2009). In between the outer rainbands of Hurricane Isabel, the 
potential temperature mixed layer depth was 400 m, while the momentum and humidity 
fluxes decreased to zero at 700 m, and the inflow layer extended to ~1 km. A similar 
depth of the inflow layer at the eyewall was reported in Montgomery et al. (2006), but the 
virtual potential temperature mixed layer was only 150 m deep. Smith et al. 2008 have 
suggested that the PBL has different characteristics in the main updraft region versus the 
inflow region outside the eyewall. In the eyewall, turbulent eddies may be transported 
upward where flow erupts abruptly out of the PBL. Under such circumstances, one might 
expect these turbulent contributions to increase the magnitude of the unresolved turbulent 
fluxes at the top of the control volume. Examining a range of relevant altitudes is a 
practical approach to addressing this uncertainty. Sensitivity tests described in Chapter 
III.E suggest that volume tops below and near the top of the inflow layer are the most 
appropriate for the budget calculation. 
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Figure 16. Root-mean square magnitudes on a logarithmie seale of the (a) momentum 
budget terms and (b) energy budget terms from the RE87 simulations with 
different flux configurations (see inset). The colors of the budget terms along 
the bottom indicate which terms may be calculated or estimated from the 
CBLAST dataset - green are known terms, red are unknown terms, yellow is 
an indirectly estimated term, and blue is the desired surface flux term. 


2. Sea-Surface Temperature 

The sea-surface temperature plays an important role in determining the enthalpy 
exchange coefficient by determining the saturation enthalpy at the sea surface 
(k* = + Lq*) and therefore must be measured accurately. Direct measurements of 

the SST were made during the Hurricane Fabian missions by AXBTs released by the 
NOAA aircraft. A comparison of the AXBT temperatures and TRMM microwave imager 
satellite estimates from the 3-day period prior to the first mission is illustrated in Figure 
17 and summarized in Table 2. The TMI-derived sea-surface temperatures are near 28.5° 
along the track throughout the analysis domain, but the AXBT temperatures are 
consistently 1-2 degrees lower at the same locations. Since decreases of the SST are 
known to occur during and after the storm passage, the splash locations of the 
measurements were plotted in a storm-relative coordinate system (Figure 18). The AXBT 
measurements were near the eyewall (red annulus) with the exception of AXBTs #4 and 
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#8. On 2 September, the two AXBTs on the right side of the eyewall are ~1 C lower than 
the one released on the left side. Given the relatively good azimuthal sampling around the 
eyewall, an average value of 27° C was chosen for this mission. On 3 September, all of 
the AXBTs near the eyewall agree to within 0.4° C, and are consistently -1.2° C lower 
than the TMI temperatures. A slightly higher value of 27.5° C was utilized on this day. It 
is interesting that the general location of Fabian on 3 September coincides with that of 
Isabel on 13 September. Since the SST is estimated at 27.5° C on both of these days, it 
appears that some of the ocean cooling that occurred during Fabian’s passage recovered 
in the 10 days between the two storms. Only a single AXBT was available on 4 
September, which was adjusted upward slightly to 28° C, given the warmer SST field 
derived from TMI in this region. 

A comprehensive analysis of the SST for the potential intensity estimates for 
Hurricane Isabel in BM08 used satellite, buoy, and airborne radiometer measurements 
(Figure 19). Comparisons of the 3-day TMI product with the AVHHR estimates are 
good, and lend confidence to the prior estimates. Given the high bias of the TMI 
estimates in the Fabian case, the 29° C on 14 September used in BM08 was adjusted 
downward slightly here to 28.5° C. The 28.5° and 27.5° C sea surface temperatures 
obtained previously for 12 and 13 September, respectively, were not modified. It is 
possible that these estimates are a little high for these calculations, given the lower 
temperatures found underneath Fabian’s eyewall, but without additional evidence they 
are the best estimates available of the SSTs. 

SST errors do not affect the retrieved enthalpy flux, but do affect the 
corresponding exchange coefficient via the surface saturation enthalpy. The percentage 
errors in Ck for 0.5° C SST error increments from both the WRF and RE87 output are 
shown in Figure 20. The two models agree quite well, which is largely expected since the 
error is not dependent on storm structure but only on the prescribed saturation enthalpy 
value. It is reasonable to assume the error is not much larger than 1° C, since all of the 
SST measurements are within a few degrees of each other, and are directly measured in 

the Fabian case. This would suggest an error of no more than 20% in the derived Ck. 
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Table 2. 


Table of SST Observations for Hurrieane Fabian 


AXBT# 

D/ HH;MM;SS 

AXBT SST 

TMI SST 

1 

2/ 17:25:24 

26.64 

28.95 

2 

2/ 20:02:29 

27.61 

28.65 

3 

2/ 20:17:40 

26.88 

28.5 

4 

2/ 20:49:46 

26.72 

28.8 

5 

3/ 17:12:06 

27.73 

28.95 

6 

3/ 17:38:25 

27.4 

28.65 

7 

3/ 19:15:23 

27.34 

28.65 

8 

3/ 20:34:07 

27.64 

28.5 

9 

4/ 19:44:36 

27.49 

29.1 


AXBT - TMI 

AXBT/TMI 

Analysis 

Difference 

Average 

SST 

-2.31 

27.795 


-1.04 

28.13 


-1.62 

27.69 


-2.08 

27.76 

27.0 

-1.22 

28.34 


-1.25 

28.025 


-1.31 

27.995 


-0.86 

28.07 

27.5 

-1.61 

28.295 

28.0 


TMI Sea Surface Temperature 3-day average ending 09/02/2003 
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Figure 17. SST derived from TRMM Mierowave Imager satellite (average SST over 31 
August to 2 September in eolor), and AXBT data released into Hurrieane 
Fabian. Traek of Hurrieane Fabian (dashed best traek, from 2 to 5 September) 
is shown for reference. 
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Figure 18. 
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Storm-relative AXBT splash locations. AXBT numbers correspond to Table 
2. Red annulus corresponds to approximate eyewall location at 30 km radius 
from the storm center. 
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Figure 19. SST derived from AVHHR satellite (average SST over 4 to 10 September in 
color), and NOAA WP-3D downward-pointing radiometer (thin line, from 
~18Z 13 and 14 September). Tracks of Hurricanes Fabian (dashed best track, 
from 2 to 5 September) and Isabel (dashed best track, with thick white, solid 
line indicating analysis periods from 16-23Z on 12 to 14 September) are 
shown for reference. (From Bell and Montgomery 2008) 
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Sensitivity to SST Errors 



Figure 20. Retrieved Ck sensitivity to errors in sea-surface temperature using RE87 and 
WRF model runs. Units are in percent for the Ck error. 

3. Circulation Centers 

It is well known that an inaccurate center can lead to errors in the Fourier 
decomposition of a storm’s kinematic and thermodynamic structure in cylindrical 
coordinates (Marks et al. 1992; Lee and Marks 2000; Bell 2006). While a transformation 
of the data into cylindrical coordinates does not change a field (except for the singular 
point at the origin), the incomplete representation of the field by a finite series that 
truncates higher-order harmonics does introduce error. This limitation is particularly 
important when only the wavenumber zero mode is retained, as aliasing to higher 
wavenumbers is manifest as errors in the derived axisymmetric winds, energy, and mass 
fields. For scalar fields, center errors are manifest as radial displacements of the 
measured quantity, such that sufficient averaging can largely eliminate random errors. 
For vector fields, a center displacement in the radial direction has the same effect as 
merely moving the observation inward or outward from the displaced center. Sensitivity 
tests indicated that this type of error leads to minimal errors in the budget retrieval. 
However, displacements of the center in the transverse direction lead to a re-partitioning 
of the winds into the tangential and radial component. 


53 




It can be shown from geometrical considerations in a polar coordinate system that 
the re-projected wind components are given by u = ur! r + ve! r and v = vr I r + ue I r, 
where e denotes the transverse error in the center, u and v are the actual radial and 
tangential wind components, r denotes the radius of the observation, and tildes indicate 

erroneous quantities with r = . Since r I r ~l for typical center errors away 

from the eye, the percentage errors are approximately given by the second terms in each 
equation. Given that v » w, it is clear that the largest errors would be found in the radial 
component on the inner side of a budget control volume. For example, a 5 km 
displacement of a 60 m s'^ tangential wind at 20 km radius would yield a 15 m s'^ error in 
radial wind. Weighting this error by mass, radius, and energy or angular momentum and 
integrating over the column/ring yields a large error in the resolved flux across the 
boundary of the control volume. In cylindrical coordinates, integration errors are 
compounded as the depth of the volume increases, which lends some support for the use 
of shallower control volumes that still exceed the depth of the boundary layer. The fact 
that the error is radially dependent means also that an artificial gradient of the radial wind 
is introduced, which produces an erroneous residual flux into the volume and inaccurate 
radial divergence. Since changes in the divergence affect vertical motion, the error also 
would be spread to the flux across the top of the control volume. All of these error 
implications add sensitivity to the choice of the control volume. 

Sensitivity tests with the RE87 model confirm the large errors described above 
when applied to an individual model output, whether instantaneous or time-averaged. 
Standard deviations of the percentage error for both drag and enthalpy coefficients were 
-25% per km center error, such that center errors of 4 km had a ±200% variation at the 
95% confidence level. It is partly for this reason that individual radial penetrations and 
dropsonde sequences were not used to calculate the budgets, since a center error can be 
significant on a single pass. Fortunately, if the center displacements are predominately 
random, then averaging of multiple passes with independently derived centers reduces the 
errors in the winds and thermodynamics, which suggests that the analytic results are 
really a worst-case scenario. However, this also indicates that reducing center errors is 
critical to obtaining meaningful results from the budget. 
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Since a variety of ways exist to define the center of a tropical cyclone, there is no 
unique tropical cyclone center, but there may be an optimal center to minimize aliasing 
associated with a truncated Fourier representation. The simplest definition of a TC center 
is where the earth-relative wind goes to zero. For a moving vortex, the wind center then 
depends on the translation speed and is not ideal for most quantitative applications. The 
most commonly used center in TC studies is the 6-hourly “best track” center published by 
the National Hurricane Center. These centers are derived from a variety of sources and 
are accurate to approximately 10 km for strong hurricanes. While some of these centers 
can be highly accurate when aircraft reconnaissance is available near the 6-hourly 
interval, the low temporal resolution of the fixes smoothes out any high spatial resolution 
movement such as trochoidal oscillations (Nolan et al. 2001). Although the best track is 
sufficiently accurate for comparison with forecast tracks, which typically are concerned 
with errors >100 km, it is not optimal for deducing axisymmetric structures. In previous 
studies using the current dataset (Montgomery et al. 2006; Bell 2006; BM08), the centers 
were determined by the method described by Willoughby and Chelmow (1982), in which 
flight-level pressure and wind observations are used to find the streamfunction minimum 
for a axisymmetric vortex in gradient wind balance. In the absence of multiple local 
minima that may arise from transient mesovortices, this dynamic center is accurate to 3 
km and is the preferred method for high-resolution centers when only flight-level data are 
available. A third center-finding method proposed by Marks et al. (1992) used Doppler 
radar data to determine the circulation center that maximizes the axisymmetric tangential 
wind at the RMW. Since this maximum is an inflection point in the radial gradient of 
tangential wind, it also maximizes the vorticity and therefore the circulation inside the 
RMW. An alternate center definition using both the gradients of angular momentum and 
energy also was considered, but was difficult to implement given the need for accurate 
gradients of thermodynamic information. In the case of a stationary, circular, vertically 
aligned vortex in gradient wind balance, the wind, dynamic, and circulation centers 
would be identical. As a TC translates and its structure departs from axisymmetry and 
vertical alignment, these centers begin to diverge. 
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For the current study, the availability of edited, dealiased Doppler radar data made 
the circulation center an attractive choice. To determine the circulation center, a two- 
dimensional “simplex” search (Neldar and Mead 1965) on a gridded dual-Doppler wind 
field (Mohr 1988) with 1.5 km horizontal grid spacing finds the point that maximizes the 
average value in an annulus around a specified radius. A 3-km annulus was used in this 
study, and an azimuthal average spanning at least 60 degrees was enforced to prevent 
small-scale wind maxima from dominating the average. The simplex search compares the 
average tangential velocity at the three vertices of a triangle whose centroid is an initial 
guess of the TC center. The algorithm uses reflection, contraction, and expansion of the 
triangle to find the maximum tangential velocity. To ensure a global maximum was 
found, the simplex search was initialized at 16 locations over 15 radii bracketing the 
estimated center and RMW, respectively. A mean center was deduced from the different 
simplex solutions at 1-km altitude for each radial penetration of the airborne radar. 
Cubic-spline interpolation of the individual centers over time was then used to create a 1- 
second track. The track was reviewed subjectively for outliers, and some minor 
adjustments were made after identifying centers that may have been biased by local wind 
maxima. The accuracy of individually derived circulation centers is ~3 km due to errors 
in the gridded dual Doppler analysis and spatial averaging used in the simplex search. 
However, the use of a broader swath of radar-derived wind data helps to eliminate centers 
associated with mesoscale vortices that distort the center fixes from the Willoughby and 
Chelmow (1982) method. Most importantly, since these centers maximize the symmetric 
tangential wind by definition, they also minimize any artificial high-order harmonics in 
the Fourier decomposition. 

An example of the center-finding method for a radial penetration of Hurricane 
Isabel at 1650 UTC 12 September 2003 is shown in Figure 21. A composite radar 
reflectivity from the lower fuselage 5 cm, C-band radar (color) reveals the nearly 
symmetric eyewall with values exceeding 40 dBZ that is superimposed on a swath of 
winds at 1-km altitude from the tail Doppler radar on the southwest side of the vortex. 
The white typhoon symbol indicates the location of the radar-derived circulation center, 
and the white circle indicates the estimated RMW. Note that the vectors are nearly 
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tangential to the white RMW circle, and that the RMW circle also follows the curvature 
of the outer edge of the eyewall reflectivity well in this region. Discrepancies between the 
circle and reflectivity elsewhere are likely due to temporal and spatial averaging of the 
higher beamwidth (e.g., coarser) lower fuselage data. The NHC best-track center 
interpolated to this time (black closed), and dynamic center from the flight-level data 
(black open) are also shown for comparison. The best track center in this case is very 
close to the circulation center, with the dynamic center displaced to the southwest in the 
low reflectivity region. It is possible that the dynamic center was influenced by one of the 
large mesovortices present in the eye at that time (Kossin and Schubert 2004; BM08). 
This comparison indicates that the radar-derived circulation center is the optimal one for 
the axisymmetric budget calculations. 

The NHC best-track and the high-resolution tracks derived from the simplex 
searches are shown for all six missions in Figure 22. As expected, the tracks are nearly 
coincident, but the higher resolution motion is apparent in the colored tracks. The largest 
discrepancy appears to be on 13 September, where the differences are likely due to the 
different center definitions, but are still well within the large eye and RMW on this day. It 
is reasonable to assume that no systematic biases are present in the tracks, such that 
random errors in the center would be largely minimized by the variational analysis 
technique. Sensitivity tests with minor variations to the observed tracks were consistent 
with this hypothesis, and did not indicate any significant changes in the retrieved 
exchange coefficients for displacements on the order of 1-3 kilometers. The error 
remaining in the exchange coefficients after averaging the observations from multiple 
radial penetrations is estimated at ±20%. 
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Figure 21. 


Isabel 09/12 16;50-16:58 UTC 
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12 16 20 24 28 32 36 40 44 


Radar analysis of Hurricane Isabel from 16:50-16:58 UTC 12 September, 
showing reflectivity composite from lower fuselage radar in color with dual 
Doppler horizontal winds at l-km altitude (vectors). Hurricane symbols 
indicate centers derived from radar (white), NHC best track (closed black), 
and aircraft winds (open black). White circle indicates radius of maximum 
wind from radar-derived center. 
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Figure 22. NHC best tracks (black and gray) and radar-derived tracks (color, see inset) 
used in this study for Hurricanes Fabian and Isabel. 

4. Gridded Analysis 

To integrate the kinematic and thermodynamic integrals numerically in Equations 
(3.15) and (3.33), the data must be available at specific grid points. For the previous TC 
momentum and energy budget studies, this was done by a manual map analysis that was 
then interpolated to regular intervals. This analysis was largely superseded by the 
development of objective analysis techniques that created a gridded data field by 
weighting the observations by their distance from a particular gridpoint. Modern data 
assimilation techniques have continued to improve the analysis of data, by including 
observational error, a priori background estimates of the atmospheric state and, in some 
cases, time-dependent probabilistic information on the background errors. 

For the current study, a variational analysis technique was developed based on the 
work of Ooyama (1987) and modem data assimilation methods. The objective of this 
method is a maximum likelihood estimate of the atmospheric state for a given set of 
observations and error estimates. A detailed derivation and description of the SAMURAI 
(Spline Analysis at Mesoscale Utilizing Radar and Aircraft Instrumentation) is given in 
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Chapter IV. Although the statistical correlation between the resulting SAMURAI 
analysis and the CBLAST observations is very good and yields high confidence in the 
derived TC structures used in the momentum and energy budgets, errors arising from the 
gridded representation of the observations must still be considered. 

Although variational analysis is superior to objective or manual analysis in a 
probabilistic sense, shortcomings still exist in the practical implementation due to 
incomplete knowledge of the observational and background errors. In a variational 
analysis, the derived structure conforms to new observations where available, and relaxes 
back to the a priori background state where no new information is provided. Given the 
relative lack of gridded, axisymmetric analyses of mature TC structure, the error 
estimates for the background state are uncertain. For the current study, a Gaussian 
background error covariance was assumed with an influence length scale based on the 
observation density. The background errors were purposefully set high in order for the 
analysis to conform primarily to the observations where available, but this has the 
detrimental side effect of making the analysis unconstrained in data-poor regions. 
Similarly, the background error length scale acts as both an effective distance for the 
observations and as a spatial filter; a large value helps spread the information provided by 
an observation but also tends to overly smooth the analysis. An appropriate balance must 
therefore be struck between these trade-offs. This balance was assessed by varying the 
specified length scale and subjectively examining the resulting gridded fields of the 
numerical model output and the CBLAST analyses. 

With the exception of the transverse streamfunction y/, the majority of the gridded 
fields were found to be relatively insensitive to the specific choice of the error length 
scale other than the level of detail. Since the secondary circulation depends on the 
derivatives of this streamfunction, even relatively smooth fields of y/ can produce 
substantial variations in u and w in regions where the derivative is not properly 
constrained by either the observations or the background field. This deficiency was found 
to be most significant in the lowest 100 meters of the analysis domain, where the signal 
from many dropsondes was lost and no Doppler radar information was available. An 

example of the sensitivity of the radial wind at the surface to dropsonde coverage and 
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three background error length scales is shown in Figure 23. In the main eyewall region, 
the radial wind profiles are very similar and provide a good least-squares fit to the 
observations. In the data gap regions near the eye (10-20 km radius) and outside the 
eyewall (>50 km radius), the solutions diverge significantly. The most radical is the 4x 
(e.g., 400 m vertical, 4 km radial) length scale, which has a large oscillation in the radial 
wind in the data-poor region. However, this was eliminated largely by using the 6x filter. 
The differences between the 4, 6, and 8x radial winds are pronounced at the outer 
boundary, with nearly a 40 m s’' difference where no observations are available to 
constrain the solution. This is partially because the absolute errors in Xj/ increase with 
radius as the magnitude of the streamfunction increases, and partially due to the lack of 
an ideal vertical boundary condition on the streamfunction near the surface. Since y/ = 0 
at the surface constrains the value but not the derivative, setting the first or second 
derivative to zero introduces an artificial structure in the radial wind. In this study, a 
third-order derivative constraint is used to damp oscillations of y/ near the surface, but the 
magnitude of this constraint is constant over the domain. These tests indicated that the 6x 
background error length scale is a good trade-off between too much smoothing and the 
data density constraints, and is adequate for the current study. More research on the 
proper boundary condition for the cubic spline analysis of the streamfunction is 
recommended for future applications, with alternative spline coefficients, a variable 
derivative constraint, or anisotropic background error covariance possibly improving the 
near-surface solutions in data-poor regions. 

An additional factor that improved the analysis was the incorporation of 
Doppler radar data. In BM08, the Barnes objective analysis of Hurricane Isabel used 
the equivalent of a ~10x filter to damp noise and spread the dropsonde information 
across data gaps. By incorporating radar data, the dramatic increase in both azimuthal 
and radial data density for the kinematic variables adds significant value to the analysis 
and helps reduce the amount of filtering required. Without the radar data, the surface 
radial wind sensitivity (Figure 23) extends through a much deeper layer. Ultimately, 
higher quality and quantity of observations lead to better analysis, such that the highest 
confidence in the gridded fields was in the eyewall region where the best data coverage 
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was available. The budget control volumes were therefore centered on the regions with 
both the best observations and the highest wind speeds. 

As an additional test of the errors associated with the gridded analysis, the RE87 
model output was treated as simulated observations and analyzed with the identical 
SAMURAI configuration used for the real observations. A comparison of the “midO.7” 
wind field at 168 hours is shown in Figure 24. The analyzed and the raw model fields 
compare well, and the structure is very similar to that in the “low” output at the same 
simulation time (c.f.. Figure 11). Whereas the raw model and analyzed tangential winds 
are nearly identical, the analyzed radial winds are noticeably smoother than the raw 
model output, especially in the peak magnitude of the inflow and outflow at the RMW. 
The effect of the boundaries is also apparent, with a broader region of strong outflow at 4 
km in the analyzed field and a bend in the radial wind contours near the surface. Since the 
lowest model was 125 m, the specification of winds near the surface is arbitrary in this 
case and is not necessarily representative of an error. In general, the analysis faithfully 
reproduces the kinematic structure with a 0.9984 statistical correlation, a ~1% RMS 
difference of the variables other than the radial and vertical wind. There is a -28% RMS 
difference for the radial and vertical wind, with most of the differences occurring near the 
boundaries. The 6x filter length is perhaps too heavy for a “perfect” data distribution, but 
is a fair trade-off with the CBFAST data. Given the uncertainties described above, it is 
estimated that errors arising from the gridded analysis procedure contribute less than a 
20% error in the derived exchange coefficients. 
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Surface radial wind filter sensitivity 
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Figure 23. Surface radial wind sensitivity to background error length scale and gaps in 
the dropsonde data on 13 September. The dropsonde wind observations 
below 50-m altitude on 13 September are marked by squares, and three 
analysis curves with different background error covariance length scales are 
shown in color (see inset). 
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(a) Original RE87'mid0.7'output at 168 hours 
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Figure 24. RE87 raw model output from “midO.7” simulation at 168 hours (panel a) 
eompared to SAMURAI analyzed fields (panel b). Tangential winds (m s'^ 
shaded, seale at bottom), and radial winds (5 m s'^ eontours) with dashed 
contours indicating inflow and solid contours indicating outflow. The thick 
black contour indicates zero radial wind. 

E. DERIVED EXCHANGE COEFFICIENT SENSITIVITY USING 

SIMULATED DATA 

The previous analysis of the (i) unresolved budget term magnitudes and (ii) 
SAMURAI gridded analysis method suggest that errors in the derived bulk exchange 
coefficients could be significant even if the SST and circulation centers were known 
exactly. This motivates the following test of the methodology using simulated 
observations from the numerical model output analyzed by SAMURAI with no 
instrument or center errors. Using this test, the following questions are addressed: 

1. How precisely can one estimate the exchange coefficients? 

2. Is there an optimal geometrical configuration (width and height) for the 
budget control volume? 
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To more accurately estimate the minimum uncertainty in the coefficients given 
“perfect” observations, a series of sensitivity tests was conducted to guide the 
specification of the optimal configuration of the budget control volume when real data 
will be used. The 6-hourly averaged output from 168-288 hours from all five of the RE87 
runs, and the hourly-averaged “figure-4” analysis from 144-168 hours of the WRF 
simulation were analyzed with the SAMURAI program using the same configuration as 
done for the real data. A wide variety of control volumes were then tested and compared 
with the known exchange coefficients. 

One of the first steps was to determine the sensitivity to magnitude of the 
exchange coefficient using the re-analyzed RE87 output as shown in Figure 25. The large 
quantity of model output and extensive variety of the control volume parameter space 
(e.g., variable inner radius, width, and depth) yielded -21,500 retrievals. From this test, 
generally more negative percentage errors are found for the “high” drag coefficients 
(>3.5x10'^), since only a few samples are above zero (Fig. 25a). A slight linear 
correlation exists in the largest enthalpy coefficients (Fig. 25b), which generally 
corresponds to increasing average wind speed in the control volume. For the lower 
magnitude exchange coefficients, no particular dependence is found on the magnitude of 
the coefficient, with numerous samples containing both high and low percentage errors. 
Although the percentage error in Figure 25 seems to be a useful metric for examining the 
problem, absolute errors for the higher magnitude coefficients may be on the same order 
as the lower values of the exchange coefficients. Since a systematic bias clearly exists for 
the “high” exchange coefficients, these were removed from further consideration to 
ensure the magnitudes correspond more closely to the actual retrieved magnitudes from 
the CBFAST dataset. A total of -12,900 retrievals remained after the “high” simulations 
were removed. This subset of the data does not have a systematic low bias, and has 
percentage errors that correspond to absolute errors with similar magnitudes. 

The sensitivity of the drag coefficient to the control volume geometry for the 
“low” and “mid” simulation samples is shown in Figure 26. The sensitivity to the control 
volume depth (Fig. 26a) has a parabolic trend with an apparent low bias for almost all 
depths. However, a positive peak in the distribution is found around 1500 m, which is just 
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below the inflow top for these simulations. The sensitivity to the control volume width 
has a decreasing spread of the uncertainty with wider volumes, but also has a consistently 
low bias at the widest annulus of 22 km. A useful metric for combining these two 
parameters is shown in Fig. 26c, where the width/depth “aspect ratio” defines the shape 
of the volume. The spread of the uncertainty decreases from low aspect ratios (e.g., taller, 
thinner volumes) to high (e.g., shorter, wider volumes), but also indicates a consistently 
low bias at the higher ratios. The sensitivity to the location of the inner edge of the 
control volume (normalized by the RMW) has a general low bias but with a peak around 
0.8 in the distribution. 

The RE87 sensitivity of the enthalpy exchange coefficient (Figure 27) has similar 
patterns with respect to the depth, width, and aspect ratio of the control volume (Figures 
27a-c), and to the location of the inner radius (Figure 27d). The parabolic shape of the 
depth sensitivity distribution peaks at a slightly higher height, and is more symmetric 
about the zero error baseline. The other three panels indicate a trend toward a low bias at 
the upper end of the tested ranges, but not as large of a low bias as the drag coefficient 
distributions at the same values. 

The control volume sensitivity tests for the drag coefficient (Figure 28) and the 
enthalpy exchange coefficient (Figure 29) with the WRF simulation have similar trends 
with respect to increasing depth, width, aspect ratio, and inner radius as in the RE87 
results, but with different bias characteristics. A parabolic shape is evident in the 
sensitivity to the depth of the control volumes, with peaks at a lower altitude near 1000 m 
and 400 m for Cd (Eigure 28a) and Ck (Eigure 29a), respectively. The percentage errors 
become more negative for increasing volume width (Eigures 29a and 29b). The aspect 
ratio sensitivity (Eigures 28c and 29c) appears to have the most similarity with the RE87 
results, and indicates a decreasing spread of uncertainty at higher ratios. However, a 
distinct difference exists in the overall bias of the WRE results, with most of the Cd 
sensitivity tests having positive errors and the Ck sensitivity tests having mostly negative 
errors. Since the WRE model had two sigma layers below 100 m, an additional sensitivity 
test was performed by varying the lowest level of the control volume (e.g., zi in 
Equations (3.15) and (3.33)). This test indicated a general tendency for slightly lower 
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exchange coefficients with zi equal to 100 m versus at the surface. However, no 
discernable high or low bias existed for either level (not shown). 

One of the most significant findings from this series of tests is that the errors were 
not uniform from model to model. No single control volume was found that minimized 
both the bias and spread of uncertainty for both the drag and enthalpy exchange 
coefficients in both models. It is probable that the magnitude of the individual unresolved 
budget terms is different in the 3D WRF simulations versus the 2D RE87 simulations, but 
the overall errors were comparable. To determine whether this was primarily due to the 
differences in the 2D versus 3D simulations, or some other facet of the simulations, an 
additional sensitivity test varying the vertical mixing length in the RE87 model was 
performed. This parameter partially controls the size of the turbulent eddies in the RE87 
simulations, and was set to 200 m in the five original simulations. Towering this value to 
100 m for the “low” flux configuration produced a simulation with a similar overall 
intensity, but with a slightly shallower inflow layer. The drag coefficient sensitivity to 
control volume depth is compared for these two simulations and the WRE model in 
Eigure 30. Although the two RE87 simulations have similar error profiles, a shift to a 
shallower peak in the parabolic profile is found for the 100 m vertical mixing length 
simulation. This suggests that at least some of the error differences are due to differences 
in the simulated structure, and not solely the 2D versus 3D geometry of the simulations. 

A common feature of all the simulations is a trend toward a decreased spread of 
the percentage error for larger aspect ratios. This can be interpreted as an increased 
numerical stability as a larger integrated surface flux term is solved with wider control 
volumes. A trade-off then exists since wider volumes must account for the variability of 
the wind speed with radius, but given the larger numerical sensitivity of the exchange 
coefficients retrieved from thinner volumes, these wider volumes appear to be necessary. 
The sensitivity tests also indicate a parabolic error trend for increasing depth of the 
control volume, with a peak in the distribution just below the top of the inflow layer. 
Although the peaks were not necessarily correlated with a zero bias, they are believed to 
be related to the levels at which the vertical eddy fluxes become small (i.e., the top of the 
boundary layer). These vertical eddy fluxes may not go to zero like in the rainband region 
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as studied by Zhang et al. (2009). However, it is reasonable to assume that the vertical 
eddy fluxes are smaller at or near the top of the boundary layer compared to the surface 
fluxes. Sensitivity analysis of the retrieved fluxes from the observations show similar 
trends in the magnitude of the fluxes with depth of the control volume. The parabolic 
profile, as opposed to an increase with depth followed by a leveling off, suggests an 
accumulation of errors in the integrated energy or momentum flux through the sides of 
the control volume with increasing height. Therefore, both physical and numerical 
reasons exist to suggest that control volumes with tops near the top of the 
boundary/inflow layer and a width that includes a sufficiently large integrated surface 
flux are the best choice for minimizing error. 

Since the volume depth is constrained by the inflow layer, the aspect ratio 
sensitivity gives reasonable guidance on the volume width necessary to reduce the 
uncertainty of the retrieved Ck and Cd. These results suggest a minimum aspect ratio of 
~20 before the spread of the percentage error decreases to ~±50%, which would 
correspond to a width of 20 km for a 1-km deep volume. The data distribution ultimately 
constrains the maximum width of the control volume because of the uncertainties in the 
structure arising from fewer observations outside of the eyewall region. Therefore, an 
aspect ratio range of 20-30 is a reasonable range of the control volume parameter space 
given the observed inflow depths and data distribution. 

A summary of the WRF and RE87 error mean and standard deviation filtered to 
these aspect ratios is shown in Figure 31. The absolute error (left panel) indicates the 
smallest bias is found in the retrieved WRF Cd and RE87 Ck, with standard devations of 
0.76x10' and 0.61x10' , respectively. A distinct low bias is found in both the RE87 Cd 
and WRF Ck, with standard deviations of the error of 0.49 and 1.04, respectively. In 
terms of percentage error (right panel), these standard deviations correspond to 21-55%, 
with biases ranging from +12 to -59%. The largest biases or standard deviations are not 
confined to a particular model or exchange coefficient, but appear to be relatively 
independent. Although the simulations are idealized and the structures are only semi- 
realistic, it is believed that these results can provide useful guidance on the error 
characteristics associated with a particular set of control volumes, but do not specifically 
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describe the biases of the real dataset. However, the significant bias and uncertainty of 
these results do indicate that the methodology is sensitive to the control volume 
specification, even within this range. It is difficult to say how much of this error is related 
to the models versus the methodology. It is reasonable to assume that some of this error is 
due to noise in the models, given the temporal fluctuations in tangential wind apparent in 
Figure 10. However, the non-negligible magnitude of the unresolved budget terms is 
likely the main source of error. The errors associated with the control volume dimensions 
are therefore estimated conservatively at ±50%, but it is acknowledged that they may be 
higher than this. 

A summary of the errors identified by the sensitivity tests is given in Table 3. The 
potential errors are those that could be present if no steps were taken to address the 
problem, with the estimated errors those that are still present after the mitigation. The 
largest source of error appears to be the neglect of the unresolved budget terms, 
specifically the volume-integrated tendency and vertical eddy fluxes at the top of the 
control volume. This error is estimated at ±50% but, as mentioned above, this may be an 
under-estimate. The remaining errors all are estimated to be on the order of ±20% after 
mitigation. Fortunately, these errors are independent and random, and there is no reason 
to expect that they would be cumulative. These percentages are reasonable estimates 
based on the analytical and numerical tests performed with the model data. These tests 
also have identified the areas for which additional effort was required to ensure the best 
possible results from the observational dataset. If left unconstrained, the potential errors 
listed easily overwhelm the magnitudes of the retrieved exchange coefficients. The errors 
are non-negligible even after mitigation, but still are believed to allow meaningful 
interpretation of the results. 
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Figure 25. 


RE87 Error vs. Magnitude 
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RE87 sensitivity to the magnitude of the exehange eoeffieients. Clusters 
indieate similar magnitude eoeffieients in eaeh of the five simulations. Eaeh 
dot represents the pereentage error of an individual retrieval of an exchange 
coefficient using a single control volume at a single (6-hourly or hourly) time 
interval from one of the five RE87 simulations. 
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Figure 26. 
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The percentage error in the retrieved Cd from the RE87 simulations obtained 
by varying (a) the depth of the control volume, (b) the width of the control 
volume, (c) the aspect ratio (defined as the width/depth), and (d) the location 
of the inner radius normalized by the radius of maximum winds. Dots are the 
same as in Figure 25. 
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Figure 27. 


RE87 Ck sensitivity to eontrol volume specification. Symbols and panels as 
in Figure 26, except for Ck- 
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Figure 28. WRF Cd sensitivity to eontrol volume speeifieation. Symbols and panels as in 
Figure 26, exeept for WRF Cd- 
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Figure 29. 


WRF Ck sensitivity to control volume specification. Symbols and panels as in 
Figure 26, except for WRF Ck- 
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Figure 30. Percentage error of Cd versus control volume depth for RE87 “low” 
simulation using 200 m (green) and 100 m (red) vertical mixing length and 
WRF (blue). 
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Figure 31. (a) Absolute and (b) percentage mean bias and standard deviation for control 

volume aspect ratios of 20-30 for combined WRF and RE87 simulations (see 
inset). 
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Table 3. Summary of estimated errors in budget retrieval. Error percentages are valid for Cd 

values from 1.9-2.4 xlO'^ and Ck values from 1.1-2.4 xlO'^. 


Error Source 

Potential Error 

Estimated 

Uncertainty 

Mitigation 

Mass continuity 

residual 

>500% 

0 

Use streamfunction 

for analysis 

Shear flux term 

-60% 

±20% 

Estimate term using 

derived Cd 

Other unresolved 

budget terms 

±200% 

±50% 

Constrain control 

volume size 

SST Errors 

-20% per 1° C error 

±20% 

Use AXBT data 

when available 

Center Errors 

±25% per km error for 

un-averaged analysis 

±20% 

Use high-resolution 

radar circulation 

centers and 

averaging 

Analysis Errors 

Earge with hand or 

objective analysis 

±20% 

SAMURAI with 6x 

filtering 
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IV. SAMURAI ANALYSIS TECHNIQUE 


Improvements in objective analysis techniques for gridding meteorological 
observations have been closely coupled to the increasing volume and variety of data 
sources used in NWP. The data assimilation problem, in which the best initial conditions 
for a weather forecast must be derived from sparse observations, is essentially identical to 
the analysis problem, in which the best estimate of the atmospheric structure must be 
derived to test scientific hypotheses. In both cases, the researcher or forecaster is 
attempting to determine the “true” state of the atmosphere from a limited set of 
observations that contain errors that result from sampling, instrumentation, 
representation, and inversion. An analysis that is used to produce a forecast also requires 
balance enforcement necessary to minimize oscillations of the primitive equations, versus 
a stricter adherence to the observations for purely diagnostic studies. 

Early objective analysis techniques, such as those proposed by Cressman (1959) 
and Barnes (1973), emphasized the primacy of the observations by minimizing the 
difference between a gridpoint value and nearby observations weighted by their spatial 
distance. These well-tested algorithms still are commonly employed for diagnostic 
studies since they provide a simple, functional methodology for deriving gridded 
atmospheric structure from irregularly spaced data. These methods have been largely 
superseded within the NWP community by variational and Kalman filter techniques, 
which is motivated in part by the numerical balance requirements, but is also due to 
advancements in the application of probabilistic theory to the data assimilation problem. 
For the current study, these advancements are adopted within a mesoscale analysis 
framework in order to deduce the most probable axisymmetric TC state for momentum 
and energy budget calculations. The Spline Analysis at Mesoscale Utilizing Radar and 
Aircraft Instrumentation (SAMURAI) software program has been developed specifically 
for this study, with the goal of obtaining an objective analysis that can incorporate 
aircraft data and maximize the advantages of working within the axisymmetric, 
cylindrical geometry. 
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First, a variational (VAR) technique is derived from Bayesian first principles to 
determine the maximum likelihood estimate of the kinematic and thermodynamic 
structure. While 3D (or 2D in this case) VAR approaches are not new, the specific design 
of this implementation was chosen to provide an optimal framework for mesoscale 
analysis, and especially for tropical cyclones. Given the need for high-quality vortex and 
mesoscale analysis, particularly for initialization of numerical models, the SAMURAI 
software can be extended naturally into asymmetric cylindrical and Cartesian geometries 
in the future. 


A. 


DERIVATION OF MAXIMUM LIKELIHOOD ESTIMATE 


The derivation starts from a Bayesian approach following Lorenc (1986). Bayes’ 
theorem states that the posterior probability of a state A, given that state B exists, is 
proportional to the prior probability of state A multiplied by the probability of state B if 
A exists. That is: 


P(A I B) oc P(B I A)P(A) 


(4.1) 


where state A is defined to be the correct estimate of the atmospheric state given some 
prior observation information. In this case, the analysis state vector x is equivalent to the 
true state of the atmosphere Xt. The prior probability that this estimate is true is then 

R(A) = R(x = x,) (-4 2) 


Since some error in this estimate is likely, this error is denoted by x^-x^-x, and 
the spatial covariance of these “background errors” is then B = [x^x], ]. Here, the 
fundamental assumption is that the errors are unbiased and the probability distribution is 
Gaussian about the background estimate. This assumption can be justified by the 
expectation that the errors are random, as opposed to being systematic. The Central Limit 
Theorem for a continuous, multivariate system yields: 


P{x = x,) = 


1 


G'JIk 


exp 


^ (x - pf ^ 


1 


|B|2 yfln 


-exp 


( \ ^ j A 

V 2 y 


(4.3) 
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where |B| is the determinant of the covariance matrix equivalent to the standard deviation 
(7 in the univariate distribution, and Xb is the background state estimate equivalent to the 
sample mean ji. State B is then defined to be the state where the set of observations are 
obtained by sampling the (true) atmospheric state. That is, P(B) = P(y = y, I x = x,). 
Assuming that perfect observations could be obtained by some forward model h(x) that 
maps from the true state to observational space, then the error is denoted by 
yi;-y,- h{x^) and the covariance is R = [y^yj]. These errors encompass instrument and 

representation errors. This distribution is again assumed to be unbiased and Gaussian 
which yields a similar probability distribution function: 


P{y = y^\x = X^)- 


—T - exp\ --(y^ -/i(x,))^R ’(y, -/i(xj) 

|r|2V^ ^ ^ 


Given these two probability distributions, Bayes’ rule is applied to derive the 
probability of the posterior distribution of state A given state B; e.g., the probability that 
the atmospheric state is correctly estimated given that measurements of the true state 
have been obtained. Substituting the PDFs into Bayes’ rule yields 

( 1 ^ 

P{x = x,\y = y,) = C-exp - - (x, - x^ f B”’ (x, - x J + (y, - /i(x, R”’ (y, - h(x, )) 

^ ^ ^(4.5) 


where C is a proportionality constant given by 

P(B) 

I 

rU |b|2 



(4.6) 


Taking the natural logarithm of both sides of (4.5) yields 


lnP(x = xJy = y,) = lnC- -(x,-x^) B" (x, -xJ + (y, -/?(x,)) R~iy,-h(x,)) 

V ^ 


Thus, the probability that the atmospheric state has been correctly estimated is 
maximized when the variable term on the right side of Equation (4.7) is at a minimum. 
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Defining the term in parentheses on the right side of (4.7) as a cost function J(x), and 
dropping the subscript t for clarity, the maximum likelihood estimate (or identically, the 
minimum variance) is achieved when 

min{/(v)}= min{jj(x) + /^(x)}= 

min (x - X (x - XJ + (y - /z(x))^ (y - h{x)) | 

J (4.8) 


This minimization is achieved when 

V/(x) = B-' (X - X J + H"R-’ (y - h{x)) = 0 (4 9 ^ 

where is the linearized Jacobian matrix of h(x). Note that if h(x) is linear, the 
Bayesian solution is equivalent to the least-squares approach. 

The form of the cost function in (4.8) and (4.9) is useful, but requires the storage 
and inversion of the background error covariance matrix. This is an N xN matrix, where 
N = ix j X kx n is the dimensional space of the analysis given by the product of the 

dimensions of the spatial vectors and number of variables. This inversion can be avoided 
by recasting this equation in an incremental form, which provides an improved structure 
for estimating the effects of this matrix without requiring its full storage (Courtier et al. 
1994; Huang et al. 2000). Suppose that the maximum likelihood estimate can be obtained 
by an incremental update to the initial background estimate, such that 

x = Xi^ + dx ( 4 . 10 ) 

Substituting (4.10) into the cost function (4.8): 

J(5x) = - {5xf B"’ {5x) + {h{dx) - df R-' {h{x) - d) 

2 (4.11) 


where d = y-h{Xf^) is termed the “innovation vector.” To avoid having to invert B, a 
symmetric matrix C is defined that has the same eigenvectors as B and eigenvalues of the 

square root of B, such that B = Vb Vb = C^C (Lorenc 1997). A control variable vector 
q is defined as 
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q = Q ’^x = C ’(x-Xj) 

with the analysis increment in physical space recovered by 

dx = Cq 

Substituting (4.13) into (4.11) yields 

J{q) = ^ 

in which the gradient is derived by differentiating (4.14) with respect to q: 

VJ(q) = (I + C^H^R-'HC)g - C^H^R-'d 


(4.12) 


(4.13) 


(4.14) 


(4.15) 


This form not only avoids the inversion of B, but has several other benefits as 
well. The identity matrix in the first term of the cost function prevents the smallest 
eigenvalue of the Hessian matrix from becoming less than unity. This form reduces 
the condition number of the Hessian, and suggests the conjugate gradient algorithm 
should converge faster than the form in (4.8) (Gao et al. 2004). More importantly, the 
control variable transform allows for the use of alternate internal state vectors in different 
coordinate systems. Since the maximum likelihood estimate is highly dependent on the 
specification of B, a form is desired in which the covariance is primarily determined by 
spatial correlations instead of cross-variable correlations. Thus, variables that are highly 
coupled (i.e., u and w in the tranverse vortex circulation) can be solved dependently 
through the use of a single variable (i.e., streamfunction). Likewise, coordinate or 
spectral transforms can be applied, which further serves to isolate the background errors 
by effectively diagonalizing the matrix. Lastly, since the background error covariance 
matrix is now applied in a forward manner, it can be treated as a linear operator that can 
be modeled by a simplified function, such as a Gaussian filter. This avoids the storage 
requirements of the N xN matrix. Specific choices related to C used in the SAMURAI 
technique will be described in the following section. 
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B. AXISYMMETRIC SAMURAI IMPLEMENTATION 


The preceding general discussion is now refined for use with axisymmetric 
hurricane structure and available CBLAST observations. Recall that the probability 
distributions associated with the background and observation errors are assumed to be 
Gaussian and unbiased, and the inversion of B has been avoided by use of the 
incremental form of the cost function. However, some description of the C and the h 
operators and the control variables is still required. The background error correlations are 
assumed to be Gaussian and isotropic, which allows for an efficient recursive filter 
operator that replicates the effects of this correlation (Hayden and Purser 1995). For the 
solution to be the most effective, control variables are desired in which the analysis fields 
are smooth and well approximated by this spatial error correlation, which suggests the 
use of the following control variable state vector in cylindrical space: 

qir,z)^{prv,y/,h',ql,p[Y 

where p is the total density, r is radius, v is tangential wind, y/ is the transverse 
streamfunction, h is the saturated moist static energy given by C^T + Lq^ + gz , with Cp 

the dry air heat capacity at constant pressure, T the temperature, L the latent heat of 
condensation, and qv the water vapor mixing ratio, g is gravity, pa is the dry air density, 
and primes represent departures from a static, background reference state. All of the 
relevant physical variables can be recovered from this distribution, and these quantities 
are expected to vary smoothly on the chosen grid. Note that prv is chosen instead of v for 
two reasons: first, to match the physical transform of the streamfunction into mass- 
weighted velocity components, simplifying the use of radar data by a uniform weight to 
the radial velocity; and second, to approximate the absolute angular momentum 
distribution in the inner core. Thus, the first three variables form the primary set for 
deducing the relevant physical quantities for use in the budget equations. The remaining 
two quantities are used to isolate the density from the mass-weighted wind components, 
and temperature from the moist static energy. The use of perturbation quantities allows 
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the analysis fields to avoid strong vertical gradients associated with the background 
hydrostatic pressure and water vapor distributions. 


To convert the control vector q from gridded variables on a physical grid to 
analysis increments, the background error covariance is formulated as an operator 
sequence of a cubic spline transform that is followed by a recursive Gaussian filter: 

^^c^c^slslsls,s,s, ..... 


where 5 b is a transform from the physical, Gaussian micro-mesh (hereafter referred to as 
the “mish”) to a nodal representation, 5 a is a transform from the nodal representation to 
cubic spline coefficients, and Sc is the recursive Gaussian filter and background error 
inflation. In matrix form, 

5x = S,S,S,q = DF(P + Q)-'0^ . n. 


This yields a set of spline coefficients in physical space that can be used to determine the 
atmospheric variables and their spatial derivatives anywhere in the domain. Each of these 
operators is examined in more detail below. 

One of the fundamental aspects that distinguishes the SAMURAI analysis 
technique from other variational solvers is the use of a Galerkin approach, which is 
similar to the Fourier spectral transform, but uses the cubic B-spline as a basis (Ooyama 
2002). This has the advantages that the interpolation coefficients for a function or its 
derivative at an arbitrary point on the grid can be easily calculated a priori for each 
observation or balance constraint, instead of maintaining a memory intensive gridpoint 
representation of complex interpolation coefficients and finite difference operators. Since 
the basis is cubic, it is computationally efficient and continuously differentiable to second 
order. The finite element representation allows for a flexible incorporation of boundary 
conditions, which is a distinct advantage over pure spectral techniques for non-periodic 
domains. The disadvantage of this basis is that it is not orthogonal, and thus requires a 
matrix inversion to obtain the spline coefficients, but this is a fair trade-off with its other 
desirable characteristics. 
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The spline transform (in one dimension) is given by: 

a = (P + Q) b 

(4.19) 

where a is the state vector of spline coefficients and 

meM 

(4.20a) 

^ = b^ = jjJr)U(r)dr 

(4.20b) 


(4.20c) 


(4.20d) 

0 (r) = O -—— , for m E M and r e D 
y Ar J 

(4.20e) 

Here, u(r) refers to a continuous representation of any variable. 

u{r) refers to the 

discrete given values of the variable (in this case, the control variables on the mish) and 

O refers to the cubic B-spline given by: 

'-{I-W if i>a>0 

O J 

i(2-|^|r if 2>|||>1 

o 

0 if |||>2 

(4.21) 


The one-dimensional spline and its derivatives are shown in Figure 32. A third 
derivative constraint is incorporated into the Sa transform (e.g., (P + Q)”^) to minimize 
Gibb’s oscillations and spline representation errors near the Nyquist scale (2Ax). This 
also acts as a sharp, 6* order filter when used repeatedly as in the conjugate gradient 
minimization. For a regular r,z grid, the inverse of (P + Q)”^ is precomputed via a 
Cholesky decomposition. Since the spacing of the internal control mish is arbitrary, the 
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integrations needed for b and (P + Q) ^ are done with a two-point Gaussian quadrature 
wherein the control gridpoints are chosen to be the quadrature points of the evenly spaced 
r,z grid. Note that the 5/ transform (4.20a) is used in the H operator converting from 
spline coefficients to observation space. Homogenous or inhomogeneous boundary 
conditions can be incorporated easily into the transform by a direct modification of the 
spline coefficients near the edge of the domain. For this study, the second derivatives of 
all control variables were set to zero at the boundaries except for \(/ and prv. The variable 
prv was set to zero at the axis of rotation (r = 0), and \(/ was set to zero at the axis of 
rotation and the surface iz = 0). The extension to two dimensions is straightforward, and 
more details on the cubic spline transform and the effect of the derivative constraint can 
be found in Hausman (2000) and Ooyama (2002). 

The operator combination 5^ = DF is the application of the background error 
covariance matrix, where D is the standard deviation of the background error and F is the 
recursive filter given by: 

n 

M 

n 

■'=' (4.22) 

with r and s referred to as the advancing and backing steps of the filter. It can be shown 
that this filter approximates a second-order Gaussian filter with minimal computational 
cost. This formulation, which is derived in Purser et al. (2003), improves upon previous 
multi-pass filters by increasing the order (n) of the approximation so that only a single 
forward and backward pass is required. Analytic results suggest that fourth order is 
sufficient for most applications. 

The resulting state vector after the sequence of transforms is the set of spline 
coefficients for the state vector increments on the physical radius grid. Lastly, the 
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variables need to be mapped to observational space. This is done using the two- 
dimensional form of the cubic spline: 

x(r,z)^ IE 

neN meM (4.23) 

Note that the summation occurs only over a few nodes near the physical point in practice, 
due to the finite basis representation. The basis function is replaced by the derivatives of 
0 m for the streamfunction to convert to physical variables: 


pru{r,Z)=Y, X 

neN meM 

prw{r,z) = X S 

neN meM 


(4.24) 


For aircraft and dropsonde observations, the direct measurements of 
thermodynamic scalars T, p, and relative humidity are transformed to h’, q^, and p’a in the 
initialization to simplify the minimization. Cartesian wind velocities are then transformed 
into their cylindrical projections by: 

pu = p{(u -UJx + (v- VJy) / r 
pv = p((u -UJy-(v- VJx) / r 

pw = pw (4 25) 

where tilde variables are Cartesian wind components, v and y are Cartesian distances to 
the storm center, and r is the radius to the storm center, which is moving at (Um,Vm). 
Radar data requires the projection of the analysis variables given the following 
transforms: 

pu [(jrsin(0)cos(7) -I- ycos(0)cos(7)) /r]-l- 
pv[vcos( 0 )cos( 7 ) - ysin( 0 )cos( 7 )) / r]-l- 
pw [sin( 7 )]= pV^ 

^d=K- sin( 0 ) - sin( 0 )cos( 7 ) - cos( 0 )cos( 7 ) .4 


where 0 and 7 are the radar beam azimuth and elevation, Vd is the observed Doppler 

velocity, and Wt is the estimated terminal fall speed of the precipitation derived from the 
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reflectivity. The Joss and Waldvogel (1971) reflectivity/rain relationship 
w™'" = -2.6pZ is used below the melting level, where p is a density correction 
(Beard 1985) and Z is the linear reflectivity. For this study, the top of the analysis domain 
was restricted to 2.5 km such that the rain relationship was used exclusively, but an 
ice/snow relationship (Atlas et al. 1973) also was coded for use in other applications. 

SAMURAI was coded in C++ with the cost function minimization performed by 
a conjugate gradient algorithm (Press et al. 2002). Several critical loops were 
parallelized using openMP to accelerate convergence, with identical results to serial 
execution. For the current study, the radial nodal resolution was set to 1 km, with a 
vertical resolution of 100 m. After some experimentation, a 6Ax length-scale (6 km radial 
and 600 m vertical) was chosen for the recursive filter, and a 4Ax scale length was chosen 
for the third derivative constraint. This defined the approximate spatial influence of a 
single observation as 6Ax and damped the amplitude of spectral features less than 4 times 
the grid spacing, with both operators effectively acting as diffusive, low-pass filters. 

Single observation tests using synthetic observations were performed to test and 
debug the C++ code. An example of the analysis increment using a single Doppler 
radar observation is shown in Figure 33. The instrumentation error of this test 
observation was set very low to induce a substantial analysis increment for illustrative 
purposes. Since the Doppler velocity contains information about all three components 
of motion, a Gaussian shaped response in the tangential velocity and streamfunction 
results. The analysis of a single Doppler velocity observation is underdetermined in a 
traditional objective analysis scheme, but in the variational formulation the specified 
background field provides the additional constraints needed to produce a full wind 
field. The use of the streamfunction as the control variable induces a toroidal 
circulation around the observation to conserve mass. When multiple observations are 
included, the result is the maximum likelihood estimate that satisfies the Bayesian 
probability derived from the error characteristics of the observations and background 
analysis. The background error estimates were set liberally to allow strong adherence to 
the observations and limited dependence on the background field. To ensure the results 


87 



were not overly sensitive to the specified background, a second “outer loop” 
minimization was performed using the results of the first cost function minimization as 
an updated background field. The specified background errors are shown in Table 1. 
The tangential wind, mass, and vapor mixing ratio fields from the Barnes analysis in 
BM08 were used for the Isabel analysis, with simple, analytic fields for Fabian. In all 
cases, the transverse streamfunction background field was set to zero. 

Though the error characteristics of the specific aircraft instrumentation used in 
this study are well known, the representation errors arising from the use of an 
axisymmetric coordinate system are not. Aircraft observations were first pre-processed 
into a moving, storm-relative, cylindrical coordinate system, such that the estimated 
observation errors were a combination of errors in both the instrumentation and 
representativeness in 2D-cylindrical space. A full description of the instrumentation and 
observation errors is given in Chapter V. 
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Figure 32. 
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(a) The cubie B-spline O on the normali z ed abscissa ^ and its first derivative 
O’; (b) the second and third derivatives. From Ooyama (2002) 


Single Observation Test 
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Figure 33. Single radar observation analysis increment from background state. Color 
indicates increment in tangential wind, and vectors indicate the increment in 
the secondary circulation. 
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V. RESULTS 


A. CBLAST DATASET 

Six intensive observing periods (lOPs) were eonduetion in Hurrieane Fabian from 
2-4 September, and in Hurrieane Isabel from 12-14 September, 2003 as part of the 
CBLAST and NOAA/NESDIS OCEAN WINDS experiments. National Hurrieane Center 
(NHC) best-traek intensities were estimated at greater than 120 kt during the -16-23 
UTC time period in whieh observations were eolleeted in Eabian (Eigure 34) and Isabel 
(Eigure 35), whieh makes both storms major hurrieanes on the Saffir-Simpson seale 
during the six lOPs. Two NOAA WP-3Ds (P3s), the NOAA G-IV, and United States Air 
Eoree (USAE) Cl30 aireraft eolleeted in situ flight-level and dropwindsonde 
observations, with additional Doppler radar and radiometer data obtained by the P3s only. 
The dataset is similar to that used in Montgomery et al. (2006), Bell (2006), and Bell and 
Montgomery (2008), but with the additional use of the University of Massaehusetts 
stepped frequeney mierowave radiometer (SEMR) and Doppler radar in the axisymmetrie 
eomposites, and the addition of the Eabian observations. 

Elight-level observations at one Hz resolution from the NOAA 42 aireraft used in 
this study were kindly provided by the NOAA Hurrieane Researeh Division. A 
rudimentary eorreetion for instrument wetting errors (Zipser et al. 1981; Eastin et al. 
2002) was applied to supersaturated dewpoint temperature measurements. This eorreetion 
assumes that the errors for the temperature and humidity sensors are equal in magnitude 
but opposite in sign, whieh was shown by Eastin et al. (2002) to reduce the majority of 
significant wetting errors but not completely remove the errors. 

National Center for Atmospheric Research (NCAR) Global Positioning System 
(GPS) dropwindsondes measure pressure, temperature, relative humidity (PTH), and 
horizontal wind speed at two Hz temporal resolution along a Eagrangian trajectory as the 
sondes fall at 12-15 m s'^ in the lower troposphere. This yields a vertical resolution of 
approximately 5 m. The PTH typical errors are less than 1.0 hPa, 0.2 C, and 5% 
respectively, and the wind errors are less than 2.0 m s'^ (Hock and Eranklin 1999). 
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Vertical velocity was derived by removing the estimated terminal fallspeed of the sondes 
as a function of pressure (Wang et al. 2009). All dropsondes were quality-controlled to 
remove noise and other instrument errors with the HRD Editsonde software. Bell (2006) 
provides an in-depth discussion of the dropsonde data from the CBLAST experiment, 
including issues related to humidity corrections for some of the data. 

The X-band (3.2 cm) wavelength tail Doppler radar employed the fore/aft 
scanning technique for all missions, which provided radial velocity data in a cone ~20 
from the track both fore and aft of the aircraft. The data were first corrected for 
navigational errors and manually edited to remove ocean returns, radar sidelobes, and 
other artifacts (Oye et al. 1995; Testud et al. 1995; Bosart et al. 2002). A large number of 
radial penetrations were edited for the current study, which included all of the multiple 
dropsonde sequences. The times of the edited radar data used in the analysis are shown in 
Table 4. Most of these legs were used for circulation center fixes, but some were too 
short to obtain meaningful results from the simplex search algorithm. One of the more 
challenging aspects of utilizing this radar dataset was the velocity ambiguity resulting 
from the use of a low Nyquist velocity during most of the flights. With unambiguous 
radial velocities of <20 m s‘\ the Doppler information frequently was “folded” multiple 
times in conjunction with the very strong winds in the major hurricane eyewalls. The 
Bargen and Brown (1980) algorithm was used to correctly dealias the majority of the 
velocities using the in situ aircraft flight-level wind as a reference point, but gaps and 
noise in the data made additional manual unfolding of many rays necessary. An iterative 
correction was made for poorly dealiased velocities evident in the dual-Doppler wind 
fields used for center finding. While some of the individual radar gates still may contain 
dealiasing errors due to the multiple folds in high gradient areas, they are not likely to 
affect the analysis winds significantly, given the large number of observations and 
averaging. 

The along-track resolution of the radar scans was -1.5 km, with a 75-m range 
resolution along the beam and 1.9 degree beamwidth. In an effort to: (1) provide more 
uniform spatial resolution; (2) reduce errors associated with individual gates; and (3) thin 
the voluminous dataset for the variational analysis, the radial velocities were averaged 
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along each beam over a minimum of five gates that increased in number with range as the 
beam spread. This averaging reduced the radial resolution to approximately the diameter 
of the beam, and decreased the number of observations by nearly an order of magnitude. 
The error was determined from two contributions: (1) the spectrum width of the radial 
velocity measurements, which represents the variance of the velocity within the pulse 
volume (Keeler and Ellis 2000); and (2) the estimated error in the terminal fall speed. The 
sum of these two sources with a minimum value of 2 m s‘^ was used for the error 
estimate. 

The University of Massachusetts deployed a SFMR for measuring surface wind 
speed and rain rate on the NOAA-42 aircraft. This C-band (~5 cm) wavelength, 
downward-pointing radiometer relates brightness temperatures at six different 
frequencies to surface emissivity and to empirically derived wind speeds. Since wind 
speed is a nonlinear observation operator, it was assumed here that the tangential velocity 
comprised the majority of the retrieved wind speed. If the surface inflow angle is less 
than 30 degrees, then the error caused by this assumption is no more than -12%. Given 
uncertainties in the calibration of the instrument and the simplification of the observation 
operator, the observation error was set to a relatively high value of 10 m s'\ Nevertheless, 
these SFMR wind speed observations provided an important constraint on the surface 
wind speed when near-surface dropsonde winds were not available. 

The combined instrument and representativeness errors for the observations used 
in this study are shown in Table 5. While there are likely some larger errors associated 
with individual measurements, these values are believed to be statistically accurate for 
the bulk of the observations. The analysis composites were not found to be strongly 
sensitive to the magnitude of the prescribed errors, but did exhibit a general trend toward 
smoother fields as the instrument error values were increased. If the errors were set too 
small, then the fields exhibited a tendency toward over-fitting, with a lower root mean 
square difference between the observations and analysis, but an unrealistic level of detail. 
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Figure 34. 



Date 


NOAA Tropical Prediction Center/National Hurricane Center (a) best track 
and (b) best track intensity for Hurricane Fabian. The three intensive 
observing periods on 2, 3, and 4 September are highlighted in panel (b). 
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Figure 35. 




Date 

NOAA Tropical Prediction Center/National Hurricane Center (a) best track 
and (b) best track intensity for Hurricane Isabel. Open hurricane symbol 
indicates transition to tropical storm strength, filled symbol indicates 
transition to hurricane strength, and “L” indicates extratropical transition. The 
three intensive observing periods on 12, 13, and 14 September are highlighted 
in panel (b). 
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Table 4. Analysis times for edited Doppler radar data for Fabian during 2-4 September and 

Isabel 12 - 14 September (all times UTC). 


Fabian 02 

Fabian 03 

Fabian 04 

Isabel 12 

Isabel 13 

Isabel 14 

1712-1726 

1704-1714 

1930-1949 

1650-1659 

1620-1642 

1618-1646 

1930-1941 

1839-1849 

2020-2035 

1721-1726 

1653-1709 

1716-1730 

1953-2002 

1849-1857 

2120-2139 

1828-1846 

1726-1742 

1816-1830 

2014-2027 

1929-1939 

2143-2159 

1901-1910 

1748-1807 

1849-1902 

2101-2114 

1940-1950 

2203-2219 

1943-1948 

1820-1841 

1906-1919 

2115-2123 

2010-2020 

2234-2250 

2006-2019 

1859-1913 

1936-1947 

2206-2219 

2100-2110 


2024-2036 

1916-1930 

1953-2009 

2223-2235 

2114-2124 


2038-2049 

1936-1950 

2013-2029 

2244-2252 

2204-2211 


2115-2130 

1953-2007 

2035-2050 





2026-2039 

2052-2111 





2043-2055 

2112-2125 





2102-2117 

2130-2147 






2151-2207 






2209-2225 


Table 5. Combined instrumentation and representativeness errors used in the SAMURAI 

composites. 


Instrument 

(Type) 

pu error 

pv error 

pw error 

q ’ error 

p’ error 

h ’ error 

Dropsonde 

(0) 

2 kg m'^ s'^ 

2 kg m'^ s'^ 

4 kg m'^ s'^ 

2 g kg-^ 

10 g m"^ 

5kJ 

Flight level 
(1) 

2 kg m'^ s'^ 

2 kg m'^ s'^ 

2 kg m'^ s'^ 

2 g kg-^ 

10 g/m'^ 

5kJ 

Doppler Radar 
(2) 

Min(a-i- 2*sin((|)), 2.0) kg m'^ s'^ 




SFMR 

(3) 


10 kg m'^ s'^ 
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B. DERIVED STORM STRUCTURES 


SAMURAI analyses were conducted for the six CBLAST missions using the data 
sources, error characteristics, and circulation centers described above. The radar and non¬ 
radar observation data distribution for the 2 September mission is shown in Figure 36 in 
both the radius-height and polar planes. It is clear that the radar data dominates the spatial 
distribution due to the wide swath of the tail Doppler radar. Comprehensive radar 
coverage provides good confidence in the kinematic fields throughout the domain, except 
near the eye where scatterers are limited. Excellent dropsonde coverage is also obtained 
in the main eyewall region near the RMW, with bracketing observations at flight level 
and at the surface from the SFMR. Good azimuthal sampling was achieved by the 
multiple sequences across the eyewall. Since only the flight level and dropsonde data 
provide thermodynamic data, the highest quality energy analysis is restricted to the ~20- 
50-km annulus. Note that the closed energy contours at the outer radii (> ~50 km radius) 
are likely artifacts from the lack of data in these regions, where the analysis relaxes back 
to the prescribed background state. Similarly, the thermodynamic structure inside of 20- 
km radius may be unreliable due to a data void region. Control volumes were not 
extended to these regions due to lack of confidence in the thermodynamic analysis 
outside of the main dropsonde sequences. However, the analysis was extended to the 
storm center and just beyond the main dropsonde sequence locations to avoid issues with 
the prescribed boundary conditions. 

A large number of observations was incorporated into the analysis composites. A 
total of -711,050 observations of a kinematic or thermodynamic variable were included 
in the 2 September analysis, with the largest number coming from the radar data. A bar 
chart showing the relative contributions of each instrument type for this mission is shown 
in Figure 37. Nearly 600,000 Doppler radar velocities were included, with the second 
largest contribution coming from dropsondes at -100,000. Flight level and SFMR 
observations were more limited due to their single altitude sampling. One of the 
advantages of the variational analysis is the ability to effectively combine these 
observations based on their individual error characteristics. The statistical comparison of 
the observations versus the analysis at the measurement location is shown in Figure 38. 
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The correlation is greater than 0.99, with a linear fit having only a small intercept bias 
and slope near 1. The spread of the difference between the observations and analysis is 
nearly uniform at different magnitudes. The units of the observations are m s"' for 
kinematic variables, g kg'^ for water vapor, kJ for moist static energy perturbation, and 10 

'j 

g m' for density. Since all of these variables have similar magnitudes in these units, the 
statistics are not dominated by any particular measurement. 

The prescribed background and observation errors and smoothing parameters 
control the amount of spread from the one-to-one centerline. The distribution of the 
differences is shown in Figure 38b, and illustrate an approximately normal distribution 
centered about zero, which is consistent with the uniform scatterplot in panel a. The large 
number of observations ensures that any single poor measurement will not bias the 
analysis significantly, and that the majority of the differences are within the combined 
observational error and variability of the three-dimensional observations in representing 
the axisymmetric mean. These statistics suggest a good fidelity of the analysis to the 
observations without excessive over-fitting. 

The data distributions of non-radar observations for the remainder of the missions 
are shown in Figures 39-43. Radar sampling was similar on all days to the first mission 
and is not shown. The dropsonde distributions are similar to the first mission, with the 
exception of little data outside of 40-km radius but excellent azimuthal sampling on 3 
September, and limited azimuthal sampling on 4 September. All of the Isabel missions 
are characterized by comprehensive radial and azimuthal sampling. Note that 
distributions are slightly different from those shown in BM08 due to the use of radar- 
derived circulation centers instead of flight-level winds as described in Chapter III, but 
the patterns are qualitatively the same. 

The derived wind fields for Hurricane Fabian are shown in Figure 44. The 
tangential winds depict a general weakening over the three analysis days, with a decrease 
in the depth of the strongest winds and slow weakening of the maximum tangential wind 
by the third day. The reduction in the depth and intensity of the tangential wind was 
concurrent with a reduction in the inflow depth and magnitude over the three days. The 

primary updraft also appears to have weakened during these three days. The RMW 
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remained consistently just inside 30-km radius, with a sharp gradient of tangential wind 
towards the center. Note that the low-level inflow continues past the RMW on all three 
analysis days, with the flow turning upwards near the high tangential wind gradient 
region. The absolute angular momentum and total energy are shown in Figure 45. The 
secondary circulation generally follows the contours of angular momentum above the 
inflow layer, with more tilted angular momentum surfaces outside the RMW on the 
second and third days. The energy contours also depict a toroidal shape that is similar to, 
but not congruent with, the momentum contours and secondary circulation. A reduced 
radial energy gradient is apparent on each consecutive day, which is consistent with the 
weakening trend in the kinematic variables. The oppositely directed gradients of 
momentum and energy are similar to those found in the numerical simulations in Chapter 
III.E, but the specific magnitudes and structures of the energy are notably different from 
the RE87 and WRE simulations. 

The analyzed kinematic fields for Hurricane Isabel for all three missions are 
shown in Eigure 46. Note that the inner radius shown is 15 km for clarity due to the larger 
domain, although the analysis domain was extended to the TC center. The qualitative 
features of the SAMURAI analysis are similar to the Barnes analysis used in BM08, but 
with enforced mass continuity of the secondary circulation and a notably sharper 
resolution of the boundary layer inflow. The characteristic “over-shoot” of the inflow 
past the RMW is apparent during all three missions as was the case in the Eabian 
analysis, with a sloping inflow top that coincides with the height of the maximum 
tangential wind. A gradual weakening and expansion of the tangential wind as described 
in BM08 is evident, with a well-defined outflow and updraft core just above the 
maximum tangential wind on all three days. The angular momentum and energy analyses 
are shown in Eigure 47. The consistent structure of the secondary circulation, momentum, 
and energy contours is similar to in the analyses for Eabian. The total energy structure is 
also qualitatively similar to the moist entropy structure shown by BM08 (c.f., their Eigure 
5), as would be expected by their similar thermodynamic variables. 

The ki nematic and thermodynamic structures apparent in the six analyses are 
consistent in their depiction of gradually weakening storms just past their peak intensity. 
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Differences exist in the evolution of the RMW, with Isabel expanding while Fabian remained 
nearly constant, and in their respective intensities. These six missions make an excellent 
dataset for the flux retrievals, since they are all nearly steady state, span approximately 20 m 
s’' in their peak tangential wind speeds, and have a variety of different inflow depths and 
magnitudes. 

The analysis statistics for all six missions are listed in Table 6. The number of 
observations varies from -500,000 to more than 1.7 mi llion on 14 September (due to the 
large domain required to encompass the 45 km RMW), with comparable differences between 
the analyses and observations on all days. Similar minimum and maximum differences, low 
skewness, and kurtosis values near 3 suggest nearly normal distributions. One exception to 
this statement is found on 13 September, where the high kurtosis and larger peak differences 
are most likely attributed to some faulty velocity dealiasing, but the low mean, RMS, and 
standard deviation of the differences suggest these are isolated errors. Linear correlations are 
near 0.99 for all days, with low intercept and slope values near one (not shown) for all 
distributions. These statistics suggest that the SAMURAI analyses are indeed maximum 
likelihood estimates for the TC structures in areas of sufficient data density. Since there are 
no high confidence a priori estimates of the structure available, other than the previous 
Barnes analysis with the same dataset, the data is only weakly constrained in poor data 
regions. Fortunately, the rich data regions encompass the areas required to obtain optimal 
aspect ratios of 20-30 for the control volumes. 
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Figure 36. Spatial distribution of observations on 2 September, (a) Non-radar and (b) 
radar observations in the radius-height plane, (e) Non-radar and (d) radar 
observations in the radius-theta plane. 
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Figure 37. Number of observations on 2 September for each instrument. Instrument 
types are given in Table 2. 



Figure 38. Analysis versus observations for Hurricane Fabian on 2 September, (a) 
Scatterplot of all observations and (b) histogram of differences between 
observations and analysis. Units of each observation are the same as in Table 
5. 
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(a) (b) 



Figure 39. Spatial distribution of observations on 3 September. Non-radar observations 
in the (a) radius-height plane and (b) in the radius-theta plane. 


(a) 


(b) 


Fabian 09/04 Non-radar observations 



Fabian 09/04 Non-radar observations 



Figure 40. Spatial distribution of observations on 4 September. Non-radar observations 
in the (a) radius-height plane and (b) in the radius-theta plane. 
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(a) 


(b) 


Isabel 09/12 Non-radar observations 



Figure 41. Spatial distribution of observations on 12 September. Non-radar observations 
in the (a) radius-height plane and (b) in the radius-theta plane. 
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Figure 42. 


Isabel 09/13 Non-radar observations 
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Spatial distribution of observations on 13 September. Non-radar observations 
in the (a) radius-height plane and (b) in the radius-theta plane. 
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Figure 43. Spatial distribution of observations on 14 September. Non-radar observations 
in the (a) radius-height plane and (b) in the radius-theta plane. 
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Figure 44. Fabian axisymmetric tangential wind (color), radial wind (contour), and 
secondary circulation (vectors) from 2-4 September. 
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Figure 45. 


Fabian axisymmetric angular momentum (color), total energy (350 kJ kg'^ + 
1 kJ kg’' contours), and secondary circulation (vectors) from 2-4 September. 
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Figure 46. Isabel axisymmetric tangential wind (color), radial wind (contour), and 
secondary circulation (vectors) from 12-14 September. 
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Figure 47. Isabel axisymmetric angular momentum (eolor), total energy (350 kJ kg'^ + 1 
kJ kg'^ eontours), and seeondary cireulation (vectors) from 12-14 September. 


Table 6. Statistical differences between analysis and all observations. Units of observations are 

listed in Table 5. 


Mission 

2-Sep 

3-Sep 

4-Sep 

12-Sep 

13-Sep 

14-Sep 

# Obs xlO^ 

7.1105 

8.5040 

5.4717 

9.8731 

6.5533 

17.448 

Minimum 

-42.912 

-37.362 

-29.793 

-53.507 

-67.852 

-43.795 

Maximum 

41.45 

41.632 

36.762 

41.57 

71.787 

34.927 

Mean 

-0.5821 

-0.50257 

-0.29302 

-0.1334 

-0.13643 

-0.3121 

Median 

-0.5448 

-0.4461 

-0.2477 

-0.24437 

-0.069 

-0.35765 

RMS 

5.0665 

4.8648 

5.1959 

5.0483 

6.1022 

4.7359 

Skewness 

0.01482 

0.16971 

0.19735 

0.10709 

0.16345 

0.089939 

Kurtosis 

1.5122 

2.0388 

0.60577 

3.2448 

13.503 

1.4067 

Correlation 

0.99129 

0.99192 

0.98862 

0.99317 

0.99099 

0.99342 
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c. 


AIR-SEA FLUXES 


The angular momentum and total energy budgets were resolved using the 
SAMURAI-derived storm structures shown in the previous section. A range of control 
volume depths near the top of the inflow layer from 400-1000 m was utilized, with the 
lowest level at both the surface and at 100 m height. Control volume widths from 10-22 
km in 2 km increments were then prescribed, such that the aspect ratio of the volume 
remained in the 20-30 interval suggested by the sensitivity analysis. The inner radius of 
the volume was then varied by 2-km increments around a central value of 0.8 normalized 
inner radius, which yields a total of 72 flux “samples” per mission. An example of the 
control volumes used on 12 September is shown in Figure 48. These volumes spanned 
the region in and around the eyewall, such that the average wind speed in the volume was 
always above 50 m s'^ for every mission. This spectrum of control volumes adequately 
represents a reasonable range of shapes and sizes for obtaining a mean and standard 
deviation of the budget residual at the eyewall on each day. The exchange coefficient 
sensitivity to the control volume parameters in this range was similar to that derived from 
the numerical simulations, with the exception of the WRF Ck depth sensitivity. A weak 
decreasing trend in Ck with increasing volume depth was deduced from the observations, 
whereas no distinct decrease above 400 m depth was found as in the WRF simulations. 
Results with a fixed 1 km control volume top were within 10% and 20% of the mean Cd 
and Ck values over all depths, respectively. 

The retrieved surface stress values are shown in Figure 49. Although considerable 
spread exists in the individual samples, a general agreement in the stress magnitudes is 
found with a mean value of 9.4 and standard deviation of 4.6 N m' . A linear fit of the 
data has a correlation coefficient of 0.22 and only a slight upward trend with wind speed. 
Missions on 2 and 13 September have the highest stress, but also have two of the deepest 
inflow layers. The 14 September analysis also has a deep inflow layer, but with a stress 
more similar to the other missions. As will be shown below, the 12 September retrieval 
has an anomalously high Ck/Cd ratio compared to the other missions, which suggests that 
the retrieved stress is possibly too low. It is not clear that this retrieval is an outlier, but 
more likely that the scatter is due to inherent uncertainties that were revealed in the 
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sensitivity analysis. Removing this mission from the sample improves the correlation 
coefficient to 0.53, and results in a more distinct upward trend with wind speed. Dividing 
by the density integral and taking the square root yields the friction velocity, which is 
shown in Figure 50. The scatter is less than in the stress retrieval due to the square root 
dependence, with a mean value of 2.8 m s'^ and a 0.7 m s‘^ standard deviation. These 
estimates also indicate a slight upward trend with wind speed with a correlation of 0.24, 
or 0.53 with the 12 September mission removed. 

The retrieved enthalpy flux is shown in Figure 51. These are the first known 
estimates of the enthalpy flux at wind speeds greater than 50 m s'\ and likely the highest 
estimates of heat flux recorded over the ocean surface. There is less scatter than in the 
stress retrieval, and an increasing linear correlation with wind speed of 0.81. While some 
of the individual samples clearly are too low (near zero on 3 September) and too high 
(4000 W m' on 12 September), the mean enthalpy flux for each mission and trend 
appears reasonable, with a variation from 764 W m' at 52 m s' to 2189 W m' at 
72 m s'^ The reasons for the reduced scatter and better correlation in the energy budget 
compared to the momentum budget are not apparent but, given the high variability in the 
fluxes from numerical model simulations, this result also is not too disconcerting. It 
would appear that the derived thermodynamic structure was steadier on the ~6 hour 
mission timescale than the kinematic structure, but it is impossible to validate this with 
the current dataset. 

The derived drag and enthalpy exchange coefficients are shown in Figure 52. The 
relatively low scatter of both results is very good considering the potential errors, and the 
enthalpy exchange coefficients have a better agreement than do the drag coefficients, as 

a 

would be expected from the previous retrievals. The mean drag coefficient is 2.4 xlO' 
with a standard deviation of 1.1 xlO' , which corresponds to ±46% uncertainty at the 67% 
confidence level, and a ±93% uncertainty at the 95% confidence level. The mean 
enthalpy exchange coefficient is 1.0 xlO' with a standard deviation of 0.4 xlO' , which is 
a ±40% uncertainty at the 67% confidence level, and an ±80% uncertainty at the 95% 
confidence level. These percentage uncertainties are consistent with the sensitivity 
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analysis in Chapter III that indicated a -50% standard deviation primarily associated with 
the neglected budget terms. The drag coefficient results do not indicate a significant wind 
speed dependence, with only a slightly decreasing linear correlation of only 0.11. The 
enthalpy coefficient shows a slight increase with wind speed, with a slope of 
0.03 (m s'^)"' and a higher linear correlation of 0.53. Given the uncertainties and the weak 
slopes of both fits, there is no statistically significant change in either the drag or enthalpy 
coefficient in this wind speed range. 

Histograms of the Cd and Ck samples (Figure 53) suggest that the distribution of 
the samples is similar and approximately normal in both cases. Low right skewness (0.56 
and 0.45 for Cd and Ck, respectively) with a slightly rounder peak than normal and 
negative kurtosis (-0.58 and -0.12, respectively) are found for both coefficients. These 
distributions yield good confidence that the sample mean values of the distribution are 
representative of the actual mean bulk exchange coefficients at these wind speeds, with 
the majority of the variability attributable to random errors associated with the 
uncertainties examined in Chapter III. The positive skewness suggests that there may be a 
small low bias in the mean, but this uncertainty is well within the standard deviation of 
the samples. The cumulative probability distribution (Figure 53, right) graphically 
illustrates the probability that the exchange coefficients are at or below a particular 
magnitude. The drag coefficient has a steep slope above 1 xlO'^, with -50%, 70%, and 
90% probabilities that Cd is less than 2, 3, and 4 xlO' , respectively. The maximum 
estimated drag coefficient is less than 5 xlO' with 99% probability, although this outlier 
estimate is not that useful given prior estimates’ magnitude (Chapter II.B). The slope of 
the Ck distribution is much flatter than the drag coefficient distribution. The maximum 

'j 

estimated enthalpy exchange coefficient is 2 xlO' with 99% probability, with an -60% 
probability that the value is less than 1 xlO' . These probabilities assume that the 432 
samples obtained from varying the control volume over each of the six missions 
accurately represent samples from the true population distribution of exchange coefficient 
magnitudes. Since the numerical simulation sensitivity tests indicated that the 
uncertainties associated with these control volumes may be similar to those obtained with 
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the real data, and the other error sources are also primarily random over the six missions, 
this is believed to be a fair assumption. 

Dividing the mean Ck/Cd gives an average ratio of 0.4, with a variation from 0.17 
to 1.05 by adding and subtracting a single standard deviation to each coefficient in each 
direction. The ratios of Ck/Cd calculated from the different samples and the cumulative 
probability distribution are shown in Figure 54. With the exception of 12 September, the 
retrieved ratios are in relatively good agreement. The mean value of all individually 
calculated ratios is 0.48 with a standard deviation of 0.27, which is similar to the 0.4 ratio 
and uncertainty obtained by dividing the mean Ck and Cd from all the samples. A linear 
fit shows an increasing ratio with wind speed with a correlation of 0.58, but this is 
primarily due to the large ratios above 70 m s'\ It is unclear why the 12 September 
mission appears to be an outlier, but this is probably due to unresolved errors as opposed 
to a significant change in the ratio above 70 m s‘\ It cannot be ruled out that the drag 
coefficient is reduced at these wind speeds. Since the enthalpy exchange coefficient for 
this mission agrees well with the other five, the drag coefficient derived for this day 
probably is too low, although this discrepancy could be because Hurricane Isabel may 
have undergone an eyewall replacement cycle shortly after this mission, and this cycle 
was accompanied with a change in the angular momentum structure and a larger 
integrated budget tendency term. Even with this mission included, the cumulative 
probability distribution of all Ck/Cd ratios indicates an -80% probability that the ratio is 
less than 0.75, and a 93% probability it is less than 1.0. Without the 12 September values 
these probabilities increase to 91% and 97%, respectively (not shown). 
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Figure 48. 
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Control volumes used for flux retrieval on 12 September. Tangential wind 
(color) and secondary circulation (vectors) are shown for reference. 
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Figure 49. Derived surface stress (r^g) from budget retrieval. Gray dots indicate 
individual samples from different control volumes, and large symbols 
indicate mean values from each research mission. Error bars indicate one 
standard deviation in (vertical) and average surface wind speed 
(horizontal). 
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Figure 50. 
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As in Figure 49, except derived friction velocity (u*) from budget retrievals. 
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As in Figure 49, except derived enthalpy flux from budget retrievals 
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Figure 51. 



























(a) Drag Coefficient 



(b) Enthalpy Exchange Coefficient 



Figure 52. As in Figure 49, except derived bulk exchange coefficients from budget 
retrievals. 
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Figure 53. Histogram of the percentage of samples in 0.5x10'^ bins for Cd (black) and 
Ck (gray) for all retrievals (left). Cumulative probability distributions for Cd 
( black) and Ck (gray) retrievals (right). 
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Figure 54. Ratio of Ck/Cd from budget method. Left panel has symbols as in Figure 49. 

Right panel shows cumulative probability distribution from all samples. 
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VI. SUMMARY AND CONCLUSIONS 


The magnitude of surface fluxes in major hurricanes has been very difficult to 
determine, and the resulting uncertainty remains an important aspect of both research and 
operations involving tropical cyclones. For the first time, enthalpy fluxes, the bulk 
enthalpy exchange coefficient, and the ratio of Ck/Cd have been estimated at major 
hurricane wind speeds using a budget methodology. The results presented in this study 
also augment the limited stress and drag coefficient estimates above 50 m s"'. This study 
improves upon previous studies using a budget method for estimating the surface fluxes 
in hurricanes by expanding the budgets to include total energy, by using improved data 
quality and quantity, by the use of a more sophisticated analysis technique, and by 
including a comprehensive error analysis. The six CBLAST missions into major 
Hurricanes Fabian and Isabel provided a unique opportunity to apply this budget method 
at a range of surface wind speeds from 52-72 m s'' and to obtain a good statistical 
sampling of flux retrievals. The new variational analysis scheme called SAMURAI was 
developed from first principles to determine the maximum likelihood estimate of the 
axisymmetric structure for the budgets by incorporating flight level, dropsonde, SFMR, 
and Doppler radar data. Although the derived exchange coefficients have some 
uncertainty, the sensitivity analysis using simulated data has allowed for a mitigation of 
some of the very large potential errors, and provides confidence in the quantitative 
uncertainty estimates. The main sources of error identified were: residuals in the mass 
continuity equation; unresolved volume-integrated tendency terms; unresolved vertical 
eddy flux terms at the top of control volume; sea-surface temperature errors; circulation 
center errors; and analysis errors. The “shear flux” term also was found to be important in 
the energy budget, and was included through an estimate of the drag coefficient obtained 
from the angular momentum budget. A total of 432 samples from 72 control volumes 
over the six analysis days provides a good sampling of the data that takes into account 
errors associated with the budget retrieval. 

Estimates of the drag coefficient from the current research in relation to previous 
studies are shown in Figure 55. The black symbols are taken from French et al. (2007) 

119 



that represent the eddy correlation CBLAST estimates and the laboratory estimates from 
Donelan et al. (2004). The blue symbols are taken from Vickery et al. (2009) showing the 
flux-profile estimates obtained from dropsondes. The current results are shown in green, 
along with 95% confidence intervals. Although the spread in the budget estimates from 
this study is non-trivial, the retrieved magnitudes are in general agreement with the 
previous studies that the drag coefficient does not continue to increase beyond ~30 m s'\ 
At the upper end of the 95% confidence level, the mean drag coefficient from some of the 
runs is near the extrapolated Large and Pond (1981) curve (dash-dot black curve), but the 

a 

mean drag coefficient averaged over all the samples is lower at 2.4 xlO' with a ±46% 
uncertainty at the 67% confidence level. At the edge of the probability distribution, the 
93% uncertainty at the 95% confidence level does not preclude the possibility that the 
drag coefficient is near the extrapolated value, but the cumulative distribution indicates 
that the magnitude has a 90% probability of being less than 4 xlO' . 

The individual enthalpy exchange coefficient samples are shown with previous 
estimates in Figure 56. A similar degree of scatter exists in the present results as those 
from the HEXOS and CBLAST eddy correlations, and the mean value of 1.0 xlO' with 
an ±80% uncertainty at the 95% confidence level is very similar to both the field and 
laboratory results. The good 0.81 linear correlation of the enthalpy flux with wind speed, 
and consistency of the exchange coefficients from the different missions yields good 
confidence in these results. Viewed as percentage errors, the uncertainty in the drag and 
the enthalpy exchange coefficients is similar. However, the enthalpy coefficient has a 
lower uncertainty in terms of absolute error. A linear fit of the enthalpy coefficient 
indicates a slight increase with wind speed above 50 m s"', but this increase is well within 
the uncertainty range and cannot be concluded definitively. This statistical sample 
indicates that the enthalpy coefficient has an -70% probability of being equal to or less 
than 1.2x10' , which is the approximate value determined at wind speeds greater than 
15 m s'^ These results suggest that it is probable that the enthalpy exchange coefficient is 
not dependent on wind speed in hurricane conditions. Since the current results implicitly 
include the effects of sea spray, this would also suggest that spray effects do not change 
the enthalpy exchange coefficient. One hypothesis is that the spray flux has simply 
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replaced the interfacial flux at these wind speeds, which leads to a similar overall 
enthalpy flux. However, this cannot be assessed by the current methodology. 

The results from this study have several theoretical and practical implications. 
The magnitudes of the drag and enthalpy exchange coefficients above 50 m s'^ provide 
some support of the similarity hypothesis of Emanuel (2003) that the drag and enthalpy 
coefficients would remain constant at extreme wind speeds. The retrieved Ck/Cd ratios 
from the present study are shown with previous estimates in Figure 57. These results 
indicate that the ratio is likely less than the ~1.0 estimate derived by Emanuel (2003) at 
the ~28 C temperature range, and perhaps may be as low as 0.4. The uncertainty in the 
coefficients and limited sea-surface temperature range precludes testing of the hypothesis 
that the enthalpy exchange coefficient is SST-dependent. There is an -80% probability 
from these retrievals that the ratio is less than the 0.75 ratio proposed as a threshold for 
hurricane development by Emanuel (1995). Recent theoretical and computational 
analyses suggest that this discrepancy is largely due to gradient wind imbalance in the 
boundary layer (Smith et al. 2008; Bryan and Rotunno 2009). These results are also 
consistent with recent three-dimensional numerical simulations that intensified to major 
hurricane status with Ck/Cd ratios as low as 0.1 (Montgomery et al. 2010). Although 
evidence is accumulating that potential intensity may not be as sensitive to Ck/Cd as 
originally formulated, this does not diminish the importance of having accurate estimates 
of these parameters at major hurricane wind speeds. These results suggest that the lower 
bound of the E-PI range for Hurricane Isabel presented in BM08 would be the most 
accurate (c.f., their Figure 1), which yields potential intensity estimates below 50 m s'\ 
These new estimates of Ck and Cd should help to improve potential intensity theory and 
understanding of tropical cyclone intensity change. From a modeling perspective, these 
estimates should provide additional confidence in an improved physical basis for surface 
layer schemes that do not increase the drag or enthalpy coefficient as the wind speeds 
increase above 50 m s'\ A continued examination of the physics parameterizations at 
high wind speeds in WRF and COAMPS, especially the new COAMPS-TC that is 
currently under development, will be conducted as part of future research. 
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Several additional areas beyond this study deserve further attention. With this 
proof of concept in place, fluxes in other storms with sufficient research data quality 
could be examined to improve the statistical sampling and potentially expanding the wind 
speed ranges presented in this paper. Several tropical cyclones with good reconnaissance 
and research data exist that would be candidates. In order to improve the method, the 
budget error analysis indicates that two of the largest sources of error are the unresolved 
tendency term and circulation center errors. One possible approach would be to relax the 
axisymmetric constraint, and reformulate the budgets in a Cartesian coordinate system 
centered on a particular region of the storm. This would remove truncation errors 
associated with the axisymmetric cylindrical transform, but would likely increase errors 
associated with the unresolved terms. If multiple passes could be conducted through a 
region, the tendency term could be estimated on relatively short timescales. The inclusion 
of radar data to completely document the kinematics in the control volume would be 
optimal. Since it is unclear whether this approach would reduce the uncertainty, a 
comprehensive sensitivity analysis similar to the one conducted herein is recommended. 
Lastly, the SAMURAI analysis package will continue to be developed, with a natural 
extension to three dimensions and the Cartesian domain, and with the addition of other 
data sources and analysis variables that are relevant for mesoscale studies. The usefulness 
of the software for tropical cyclone data assimilation applications with numerical weather 
models also will be investigated. 



Figure 55. Mean drag coefficients from this study (green circles) compared with 
previous studies. Black symbols adapted from French et al. (2007) and blue 
symbols adapted from Vickery et al. (2009). Red line indicates measured 
(thick) and extrapolated (thin) Large and Pond (1981) drag coefficient. 
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Figure 56. Wind speed dependence of Ck from this study (green squares) compared with 
previous studies as summarized by Haus et al. (2010). ASIST laboratory 
results (blue circles) and CBLAST (red triangles) measurements shown with 
HEXOS results (gray x’s). The mean and 95% confidence intervals are 
shown in black. 




Figure 57. Wind speed dependence of Ck/Cd from this study (green squares) compared 
with previous studies as summarized by Haus et al. (2010). ASIST laboratory 
results (blue circles) and CBLAST (red triangles) measurements shown with 
HEXOS results (gray x’s). The mean and 95% confidence intervals are 
shown in black. Purple dashed line indicates 0.75 ratio. 
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